DCE/NASA/001 7-1 
NASA CR-1 74638 


N8511456 


j~~ ' J Li 




L^v, 


'3/^ 


Xi'CvT I , 


W7 


'- 0 =? U _ vii/ _ 


I =1 r [|/ 3 \ , r 

'<=? J X? J . 


L UU 


J J w? G? J U J XL XJ vs? L 


(EASA-CE- 1746 38] PECSPiiOElC AC.d 
POWE2 PLA21 SXSIEE PEkFCEfiAMCE SCD£t C£LL 
COHPUIEH PEOGEAE Pinal Beport 'ClfJ f flD 
State Oniv.) 134 p HC AC7/E£ f*, ai> ‘ d 

'91 LGcL IGA 


G3/4< 



D) 



S65-1 1456 

Unclas 

24386 


Kalil A. Alkasab and Cheng-yi Lu 
Cleveland State University 


> ■«' r “ a=>. n 3 ^ .7 •fj'’ ^ /" 

ljjLsiuv US/© 4 .} 


Prepared for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Lewis Research Center 

Cleveland, Ohio 

Under Contract NCC 3-17 


for 




fpX 


O': 

r\lXi 'j 



\T1 


^5? 


r?' r.A'’. /R\ 
v ^7 LL 7 _ L U J _ 


>. /9=r r x n 
'*aSc/ 


' 5 =$' VS? J ~ 'oVs? J 



This report was prepared as an account of work sponsored by an agency 
of the United States Government, Neither the United States Government 
nor any agency thereof, nor any of their employees, makes any warranty, 
express or implied, or assumes any legal liability or responsibility for the 
accuracy, completeness, or usefulness of any information, apparatus, 
product, or process disclosed, or represents that its use would not 
infringe priyately owned rights. Reference herein to any specific 
commercial product, process, or service by trade name, trademark, 
manufacturer, or otherwise, does not necessarily constitute or imply its 
endorsement, recommendation, or favoring by the United States 
Government or any agency thereof. The views and opinions of authors 
expressed herein do not necessarily state or reflect those of the United 
States Government or any agency thereof. 


Printed in the'United States of America 
Available from 

National Technical Information Service 
U.S. Department of Commerce 
5285 Port Royal Road 
Springfield, VA 22161 

NTIS price codes 1 

Printed copy: A07 

Microfiche copy: A01 


TCodes are used for pricing all publications. The code is determined by 
the number of pages in the publication. Information pertaining to the 
pricing codes can be found in the current issues of the following 
publications, which are generally available in most libraries: Energy 
Research Abstracts (ERA); Government Reports Announcements and Index 
(GRA and I); Scientific and Technical Abstract Reports (STAR), and 
publication, NTIS-PR-360 available from NTIS at the above address. 



DOE/NASA/OOI 7-1 
NASA CR-1 74638 


Phosphoric Acid Fuel Power Plant 
System Performance Model 
and Computer Program 


Kalil A. Alkasab and Cheng-yi Lu 
Cleveland State University 
Cleveland, Ohio 44115 


January 1984 


Prepared for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
Lewis Research Center 
Cleveland, Ohio 44135 
Under Contract NCC 3-17 


for 

U.S. DEPARTMENT OF ENERGY 
Morgantown Energy Technology Center 
Morgantown, West Virginia 26505 
Under Interagency Agreement DE-AI21-80ET17088 




TABLE OF CONTENTS 


Page 

INTRODUCTION 1 

I. SYSTEM DESCRIPTION 2 

II. PERFORMANCE MATHEMATICAL MODEL 5 

2.1 Modeling of Fuel Processing Subsystem 5 

2.1.1 Heat Exchangers 5 

2.1.2 Shift Converters 7 

2.1.3 Reformer 9 

2.1 .3.1 Lumped Model 9 

2. 1.3. 2 Distributed Model ... 11 

2.2 Modeling of Fuel Cell Stack Subsystem 17 

2.2.1 Mass and Energy Balances 17 

2.2.2 Voltage-Current Density Characteristics 18 

2.2.3 Stack Efficiency 21 

III. PERFORMANCE COMPUTER MODEL 22 

3.1 Main Program 22 

3.2 Subroutines: BURN, CDPH, COMP, COND, CONV, DIVID, DMIX,. ... 32 

ENFU, ENRE, ENSH, EQUK, FLAME, FUCE, HEPD, HEXC, 

PDFU , PDSH, PUMP, PUP, REF, SNAE, SEPAR, SHIFT, VINEW 

3.3 Subroutine KREF 57 

3.4 Program Operation 66 

3.5 Sample Problem 72 

3.6 Discussion 79 

3.7 Further Development 79 

REFERENCES 80 

LISTING OF THE PERFROMANCE COMPUTER MODEL 81 


111 




INTRODUCTION 


1 


This report has been prepared by Cleveland State University for NASA Lewis 
Research Center to record the work done and to serve as documentation of the 
computer programs prepared under contract NCC3-17. 

Under support contract C-44219-D, energy, mass, and electrochemical 
analysis in the reformer, the shift converter, and the fuel cell module were 
combined to develop a mathematical model for the performance of the phosphoric 
acid fuel cell system which is depicted in Figure 1. 

The primary objective of the work performed under contract NCC3-17 was to 
derive the mathematical model and the associated digital computer program for 
optimizing cost and electric energy output of the phosphoric acid fuel cell 
system. To achieve this objective, all equations relating to system perform- 
ance which were derived under the previous contract, were integrated into a 
computer program that determines electric output, heat generation rate, and 
the effects on system performance of such parameters as operating pressure and 
temperature, reformer heat transfer area, and hydrogen fractional utilization. 
In addition, the mathematical and associated digital computer models were 
derived for the power processor, system components and operation costs, and 
optmization of fuel usage and cost of electric energy output. 

The present report describes just the basic performance model of the fuel 
cell system, and the computer programs written for its analyses. Other 
reports are being prepared for the cost and optimization programs, which are 
hosted by the basic performance code, and for more detailed studies of 
subsystems such as the fuel cell stack, the fuel reformer, and the heat 
exchanger network optimization. 
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A listing of the steady state performance lumped model is included at the 
end of the report. It begins on page 81. 

I. SYSTEM DESCRIPTION 

As shown in Figure 1, methane which is circulated by compressor (C) is 
preheatea by heat exchanger E-l prior to mixing it with the super heated steam 
which receives its heat by passing through heat exchanger E-9. Before entering 
the reformer, the methane steam mixture is heated via heat exchangers E-2 and 
E-3. Inside the reformer, methane is catalytical ly reformed by reaction with 
excess steam to proauce carbon monoxide, carbon dioxide, and the desired pro- 
duct, hydrogen. The effluent from the reformer is cooled by flowing through 
heat exchanger E-2 before it enters the high temperature shift converter S-l. 
The function of the high temperature shift converter is to increase the hydro- 
gen concentration and to reduce the carbon monoxide concentration of the 
reformer gas effluent. The temperature of the effluent from the shift conver- 
ter S-l is then reduced by passing through heat excangers E-l, E-9 and E-6 
Defore entering the low temperature shift converter S-2. The low temperature 
shift converter further increases the hydrogen concentration by promoting the 
shift reaction at a lower operating temperature. The effluent from the low 
temperature shift converter then enters the fuel cell containing H2, CO, CH4, 
C02 and H20. The fuel cell converts inputs of hydrogen ana oxygen to DC power, 
water and heat. Oxygen is delivered to the fuel cell by air compressor A, 
which also provides air to the reformer ourner. The spent fuel from the fuel 
cell anode goes to the burner after mixing with air supplied by compressor A. 




Flow dla^rrua of CSU designnd PA FC nynten 
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Before entering the burner, the mixture is preheated by the burner 
effluent via heat exchanger E-4. The spent fuel is then burned with whatever 
additional methane is needed to proviae the thermal energy necessary for the 
reformer reaction. 

Heat generated in the fuel cell is removed by heat exchangers E-7 and E— 10 - 
Heat from heat exchanger E-7 can then be utilized in industrial heat processing 
or space heating and cooling, while exchanger E— 10 is used to preheat the water 
supplied by liquid separator Q to provide the necessary steam needed for the 
reforming process. The effluents from the burner and fuel cell cathode will 
have their water removed and separated by condenser E-5 ana liquid separator Q 
before allowing them to be exhausted to the atmosphere. 
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II. PERFORMANCE MATHEMATICAL MODEL 

The mathematical model developed provides the basis for determining fuel 
cell voltage, current, and heat generation rate in terms of such parameters as 
flow rate, fuel composition, operating temperature, operating pressure, 
reformer heat transfer parameters, and steam-methane ratio. 

In the derivation of the mathematical model, several simplifying assump- 
tions were made. These assumptions include: one-dimensional, steady state 

flow of all gas streams, ideal gas behavior of all gas components, and a 
"lumped parameter" fuel cell stack model. 

The following subsection will consider the derivation of mass and energy 
balance equations for the gases and the description of the governing equations 
for the system output characteristics (voltage, current, and heat generation). 

2.1 Modeling of Fuel Processing Subsystem 

Production of hydrogen, which is the major function of the fuel processing 
subsystem, occurs by reaction of the fuel with steam. The major components in 
this subsystem are the reformer, the high temperature shift converter, the low 
temperature shift converter, and several heat exchangers. 

2.1.1 Heat Exchanger 

A zero capacitance sensible heat exchanger is modeled in the double-pipe 


counter mode. 
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For the counter mode, given the hot and cold side inlet temperature and 
flow rates, the effectiveness is calculated for a given fixed value of the 
overall heat transfer coefficient. The mathematical description which follows 
is covered in detail in Ref. 1. 


- T ci } 




(2-1-3) 


(2-1-4) 


where Cc: capacity rate of fluid on cold side, McCpc, J/s-K 

Ch: capacity rate of fluid on hot side, MhCpc, J/s-K 

Cmax: maximum capacity rate, J/s-K 

Cmin: minimum capacity rate, J/s-K 

Cpc: specific heat of cola side fluid, J/g-K 

Cph: specific heat of hot side fluid, J/g-K 

E: heat exchanger effectiveness 

Me: fluid mass flow rate on cola side, g/s 

Mh: fluid mass flow rate on hot side, g/s 

Qy-. total heat transfer rate across heat exchanger, J/s 
Tci: cold side inlet temperature, K 
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Tco: cold side outlet temperature, K 

Thi: hot side inlet temperature, K 

Tho: hot side outlet temperature, K 

2 

UA: overall heat transfer coefficient of exchanger, J/m -s-K 

2.1.2 Shift Converters 

The function of both types of shift converters (high temperature and low 
temperature) is to further increase the hydrogen concentration and to reduce 
the carbon monoxide concentration of the reformer gas effluent. The equation, 
CO + H^O = + C0<p(water shift reaction), dominates the material changes 

in the shift converters. The methanol input fuel does not need to pass through 
shift converters because the carbon monoxide level is low. 


In the lumped model, the water shift reaction is assumed to be at equilib- 
rium at the input temperature (isothermal operation) or the average temperature 
(adiabatic operation). The material balance is 

Pco 2 PH^ (Fco 2 + x)(FH 2 + x) 

K 2 = P co PH 2 0 = (Fco - x)(FH 2 0 - x) (2-1-5) 


where K 2 : 


equilibrium constant of shift reaction at AUT 
partial pressure of component, atm 
inlet molar flow rate of component, g-mole/s 
reacted amount rate, g-mole/s 


Equation 2-1-5 can be solved for x. Newton's methoa was used in the computer 


program. 
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The energy balance equation for the gases in the shift converter includes 
the reaction and sensible enthalpies. For adiabatic the process in the shift 
converter 


S' nj ( Ah°f ) j - S ni (Ah°f)i + S nj * (Cp)j dT 

J 298 

X- fT, 

- nc ni I (Cp ) i dT = 0 

Kb J 298 


( 2 - 1 - 6 ) 


where the subscripts PS, RS correspond to the products and reactants in the 
shift converter, respectively. Tf and Ti are the final and intial tempera- 
tures of the gases, respectively. The only unknown in the equation, Tf, is 


determined iteractively. 


The Ergun equation, which estimates pressure drop caused by the flow of 
gas through dry packings, is used to determine the pressure drop in shift 
converter and reformer. The equation is (Ref. 2): 

aP = 1878 ^ ^ ^ y + 1-75) h (2-1-7) 

e dp g c p p 

where e : void fraction in bed 

y : viscosity, Kg/m-s 

dp: effective diameter of packing particle, m 

G: superficial gas mass velocity, Kg/s-m 

h: packed height, m 

3 

p : density, Kg/m 

aP: 


pressure drop, atm 
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2.1.3 Reformer 

The key component in the fuel processing subsystem is the reformer which 
catalytical ly reforms methane (methanol or naphtha) by reaction with excess 
steam to produce carbon monoxide, carbon dioxide, and the desired product, 
hydrogen. The overall reactions are: 

C n H m + 2nH 2° = nC0 2 + ^ 2n + m/2 ) H 2 
for naphtha and methane, and 

CH 3 0H + H 2 0 = C0 2 + 3H 2 

for methanol. For simplicity, methane will be the only input fuel in the 
following discussions. 

Two reactions are assumed to be the principle reforming reactions in the 
methane-reformer, they are: 

CH^ + H 2 0 = CO + 3H 2 (demethanation reaction) 

and 

CO + H 2 0 = C0 2 + h 2 (water shift reaction). 

Reference 11 lists all of the possible reactions and discusses the minimum 
steam to carbon ratio (S/C) required to avoid carbon formation. 


2. 1.3.1 Lumped Model 

In the lumped model both of the reactions, demethanation and shift 
reaction, were assumed to be at equilibrium by utilizing the respective ADT's 
of each. The equilibrium constants were determined from the temperature. The 
equilibrium expression are 
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p p3 y 

CO2 co^ 

K l = 


ch 4 ‘ h 2 o 


Y Y 

ch 4 t h 2 o 


p p u 
c °2 H 2 
~ P P 

c CO H 2 0 


Y Y 
co 2 H 2 

~Y Y 
co T H 2 0 


(demethanation) 


(water shift) 


where Kl and K2 are the equilibrium constants of demethanation and water shift 
reaction, respectively. Expressing the mole fractions as the individual molar 
flows divided by the total molar flows yields: 


(F C0 -x+y)(F H2 +x+3y) 3 P 2 
( F CH4 _y ^ F H20 _X_y ^ F T +2y ^ 


( 2 - 1 - 8 ) 


and 

(Fco ? + x)(F +x+3y) 

K 2 = ( Fco-'X+y ) ( F H2Q -x-y ) (2 " 1_9) 

where y is the conversion amount rate in the demethanation reaction and F is 
the total inlet flow rate. Equations (2-1-8) and (2-1-9) can be solved for x 
and y. Newton's method was used in the computer program. 


The quantitites involved in the energy balance will be the sensible 
enthalpies of the gases, the reaction enthalpies of the gases, and the heat 
transferred from the combustion gases to the reformer gases, Qg_ R * The 
value of Q b _ r can be determined from 


Q b _ r = UAaT m = Hout - HIN 


where, AT m is the log mean temperature defined as 


* T m = 
m 


‘ T fc - T iR)-< T a - V 

T* - T. d 
0 fc lR 

Xn — 1 — th — 
a fR 


( 2 - 1 - 10 ) 


( 2 - 1 - 11 ) 


10 
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where, T^ is the temperature of the combustion gases after leaving the 
reformer; Tj R and T^ R are the temperatures of the reformer gases before 
entering and after leaving the reformer; A is the heat transfer area; and U 
is a modified form of a heat transfer coefficient. 

Thus, from the first law of thermodynamics and equation (2-1-11), the 
energy balance for the reformer gases can be written as, 

r Tfp 

UAaT d = m i (Ahf), - T “i (^hf)i + 21 mj (Cp)j dT 

PH J PR PR ^298 

— Z — f (Cp)i dT, ... (2-1-12) 

rR J298 

where the subscripts PR and rR stand for products and reactants in the 
reformer, respectively. 

2. 1.3. 2 Distributed Model 

Kinetical analysis was used for simulation of the performance of the 
reformer. The reformer is basically a nonadi abatic, noni sothermal catalytic 
reactor that is heated on the shell side by combustion gases from burner. 

Methane will be the only input fuel considered in this model. Figure 2 shows 

its simplified scheme. 
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Combustion 



Figure 2 Simplified Reformer Diagram 

In driving the mathematical model, the following assumptions were made: 

1. The demethanation reaction is assumed to be kinetically controlled and, 
hence, occurs at a finite rate, while the water gas shift reaction is assumed 
to be equilibrium controlled. The demethanation reaction used in this model 
is slightly modified with linear combinations of the original demethanation 
reaction and shift reaction, which results in 

CH4 + 2H20 = C02 + 4H2 (2-1-13) 

In the equilibrium calculations, the demethanation reaction choice causes no 
changes in the final results. However, the kinetic consideration will cause 
the final results to vary slightly with the reaction choice. 
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2. Axial dispersion and radial gradient are negligible - plug flow condi- 
tion. Generally, if the ratio of the length of the reactor to the catalyst's 
diameter is greater than 100, the axial dispersion effect is negligible. 

3. A uniform temperature exists throughout each catalyst particle, and 
this temperature is the same as the gas temperature in that section of cata- 
lyst bed. 

4. The kinetic expression represents a global rate, and, therefore, 
neglects reactivity differences found between the inside ana outside of the 
catalyst particles. 

5. Entrance effects are negligible. 

6. Heat transfer by radiation is negligible. 

7. Since tubular reactors inside a furnace are useo commercially, it will 
be assumed that distribution of the gas to various parallel tubes is uniform 
and, hence, a single tube is sufficient for the purpose of theoretical investi- 
gations. 


8. Ideal gas behavior is assumed. 


9. The outside shell wall is adiabatic. 
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A more detailed discussion of assumptions 3 and 4 is provided in Ref. 11 by 
examination of the "internal" and "external" effectiveness factors of commer- 
cial catalysts used in the reformer. 


Mass Balance : From the generalized continuity and the assumptions, the 

kinetic mass balance is 


where V: 
c: 

ra' : 
ee : 


v dc _ _ Y a eg 

'37 e 

average velocity of fluid through the bed, m/s 

3 

g-mole of CH^ per m fluid 

raction rate, g-mole of CH^ / s— kg catalyst 

3 

density of catalyst, kg/m bed 


(2-1-14) 


Various kinetic expressions for the reforming of methane with steam have been 
proposed which could provide the rate equation (Refs. 3, 4, and 5). The 
simplest form among the proposed expressions is the first order rate expres- 
sion, which is 

-ra' = Ko e _EA/RT P CH4 (2-1-15) 

in Arrhenius form. 


where K Q : 
EA: 
R: 
T: 


Arrhenius frequency factor, g-mole/s-kg cat - atm 
activation energy, J/g-mole 
gas constant 
temperature, K 
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Unfortunately, little agreement can be found for the values of the kinetic 
parameters, some values may be three orders of magnitude different from 
others. The data from Ref. 5, using a commercial catalyst (Gindler G-56B), is 
used i n this model . 

The water gas shift reaction is assumed to be at equilibrium. The 
conversion quantity is based upon the carbon dioxide mass balance. Thus, when 
coupled with the demethanati on reaction, the water gas shift reaction proceeds 
in reverse; therefore, the shift conversion is always negative. Using these 
two reaction schemes, all of the molar flows anywhere in the reformer can be 
written in terms of the feed quantities and the conversions of the two 
reactions. 

Energy Balance : Two energy balances are required for the system: one for the 

reformer gases and one for the combustion gases. The reformer gas balance 
includes its own sensible heat change, reaction enthalpies, and heat transfer 
from the hotter combustion gases. The combustion gas balance involves 
sensible heat change and heat transfer. This translates quantitatively into 
equations (2-1-16) and (2-1-17) 

pAiV Cp ^ = (-aH : ) & + (-aH 2 ) ^ + hixdi (Tw-t) (2-1-16) 

po Vo Ao Cpo ^ = ho*ao (T-Tw) (2-1-17) 

where aH^: demethanation reaction enthalpy, J/g-mole CH^ 

aH^: water shift reaction enthalpy, J/g-mole CO 

2 

Ai : inner tube cross area, m 
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2 

hi: wall heat transfer coefficient of tube side, J/s-m -K 

Tw: wall temperature, K 

di: inner tube diameter, m 

T: combustion gas temperature, K 

do: outer tube diameter, m 

subscript o refers to the combustion gas side 

There is greater uncertainty in estimating the heat transfer coefficient at the 
wall of tube than the rate expression. The scatter in experimental data is 
very high (Refs. 2, 3, and 4). The situation will be even more complicated by 
considering the unequal stoichimetric reaction (Ref. 6). Due to Beek's recom- 
mendation (Ref. 7), the modified Thoenes-Kramers (Ref. 8) correlation should 
be used for sphere-like particles near the wall, which are used in the model: 

hi ( dp / k f ) = 2.58(Re) 1/3 (Pr) 1/3 + .094 (Re) a8 (Pr) 0,4 (2-1-18) 

where dp: equivalent particle diameter, m 

k^: thermal conductivity, J/s-m-K 

Pr: Prandtl number 

Re: partical Reynolds number 

Differential equations, (2-1-5), (2-1-14), (2-1-16), and (2-1-17), were solvea 
simultaneously with the inlet conditions as the boundary conditions. The Ergun 
equation (2-1-7) is used to evaluate the pressure drop. 
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2.2 Modeling of Fuel Cell Stack Subsystem 

In the fuel cell power section, air, in excess of the stoichiometric 
mixture, enters the cathode side of the cell, and effluents from the low 
temperature shift converter energy at the anode. The anode input contains 
CH4, H20, H2, CO and C02. In this analysis, it is assumed that a fixed 
percentage of hydrogen is consumed at the anode, and the H20 being formed 
exits the fuel cell, with the depletea air, through the cathode exit. The 
overall reaction in the fuel cell power section is 

H 2 + 1/2 0 2 = H 2 0 (2-2-1) 


2.2.1 Mass and Energy Balances 

The lumped model provides a rapid (in terms of computation time) means of 
calculating the fuel cell module output characteristics (voltage, current, and 
heat generation rate) in terms of the inputs from the fuel processing subsystem 
and the gross fuel cell design parameters such as catalyst loading. 


The mass 


where NX: 

NI: 

Imean: 

A: 

n: 

7 : 


balances of hydrogen, oxygen and water are as follows: 


NX H2 = NI h2 - (Imean A)/(n7) 

(2-2-2) 

NXq 2 = NI Q2 - (Imean A)/(2n7) 

(2-2-3) 

NX H20 = NI H20 + ( Imean A )/( n $ 

(2-2-4) 


exit flow rate of hydrogen, oxygen, or steam, g-mole/sec 

inlet flow rate of hydrogen, oxygen, or steam, g-mole/sec 

2 

mean current density. A/cm 

2 

effective area of cell plate, cm 
number of Faraday equivalents transferred 
Faraday constant 
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The energy balance for the fuel cell is 

- (Q+W e ) =211 n i ( a hf) i - 2Z {A hf) ± 

PF J J rF 

*SInJ T fF (Cp), dT-HnJ 298 (Cp) 1 dT (2 ' 2 ' 5 ’ 

PF -'298 rF ^T i? 

where the subscripts PF, rF represent the products and reactants in the fuel 
cell, respectively. TfF is the final temperature of the products and TiF is 
the initial temperature of the reactants in the fuel cell. The nj and ni are 
the species flow rates of the products and reactants, respectively. The terms 
Q and W are the rates of heat and the electrical energy generation by the fuel 
cell, respecti vely. Q is proportional to the specific heat generation Qp 
where : 

Q = N p Xn Yn Q p (2-2-6) 

and Q F = (Tpjr - V) I (2-2-7) 


where Q: 
<¥= 
V 

Xn: 

Yn: 

I: 

AHr : 


total heat generated, J/sec 

2 

heat generated per unit area of cell, J/sec cm 

number of cells 

width of cell plate, cm 

length of cell plate, cm 

2 

fuel cell current density. A/cm 
heat of reaction, J/g-mole of H ^ 


2.2.2 Voltage-Current Characteristics 

Because of the irreversibility, the voltage V for a working fuel cell is 
the difference between the open circuit voltage and the cell polarization 


terms: 
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V = E - n 

where E: Nernst potential (reversible open circuit E.M.F.) 

n: overpotential or polarization 

The reversible cell potential, E is given by the Nernst equation: 


E 


o 


E < T > * w ln 


YH 2 /PtY0 2 

yh 2 o 


with Pt: 


E o< T > 


YH 2 : 


Y0 2 : 

YH 2 0: 


total pressure, atm 

standard E.M.F. of cell at temperature T, volts 
E (T) = 1.261-0.00025 T, T, K (Ref. 9) 


mean mole fraction of hydrogen at anode 
mean mole fraction of oxygen at cathode 
mean mole fraction of water vapor at cathode 
The polarization term n consists of four components, 


( 2 - 2 - 8 ) 


(2-2-9) 


where na: 
nr: 
nd : 
nco: 


and 


n = na + nr + nd + nco (2-2-10) 

activation poparization at cathode, volts 
resistance polarization, volts 
diffusion polarization, volts 

activation polarization at anode due to co poisoning of 
catalyst, volts 


RT , i 

na " ^oZF ln ( i o) ( SA) (CL) (CU) 


( 2 - 2 - 11 ) 
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with cC. o: transfer coefficient 

2 

i: current density, mA/cm 

2 

io: exchange current density of cathode, mA/cm 

2 

SA: specific catalyst surface area, cm / g 

2 

CL: catalyst loading on cathode, g/cm 

CU: catalyst utilization factor 

The exchange current is a function of the acid concentration, temperature, and 
partial pressure of the oxygen. The acid concentration is a function of the 
water vapor partial pressure which permits correlaion of io as a function of 
Y02, YH20, and T. An empirical fit is 

io = 232.7 (PtY02) 0 ' 8 (PtYH20) 0,4377 exp (-6652/T) (2-2-12) 

The resistance polarization is 

nr = ir 

2 

where r: specific cell resistance, ohm-cm . 

The expression of nco was chosen to have strong temperature dependence, 
be directly proportional to Yco, and have a logarithmic dependence on i, iao, 
and catalyst effective area. The resulting expression (Ref. 9) is 

nco = 0.0782PtYco exp 9190 (y - In CLa SA CU iao ( 2 ~ 2 - 13 ) 
where CLa: anode catalyst loading, mg 

2 

iao: anode exchange current, mA/cm 

Diffusion polarization has been neglected here because it is significant 
only at very high current densities. 
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2.2.3 Stack Efficiency 

The efficiency of the fuel cell to convert chemical energy to electrical 


energy, ep^, can be writen as (Ref. 10): 

e FC = E V e I e TH e H’ ••• ^ 2 “ 
where the voltage efficiency e^, the current efficiency ej, the thermo- 
dynamic efficiency e TH , and the heating value efficiency e H , are defined 
as follows: 



(2-2-15) 



(2-2-16) 


iG r 

e TH = W 


e H = aH7 
c 


(2-2-17) 

(2-2-18) 


where V and I are the operating voltage and current, respectively, E is the 
fuel cell equilibrium potential, Ip is the amount of current produced by a 
reaction, AG r is Gibb's free energy change, aH^, is lower heat of combus- 
tion of fuel cell feed, and a h f is the enthalpy change at fuel cell condi- 
tions of h 2 + \ 0 2 » H^O. 
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III. PERFORMANCE COMPUTER MODEL 

Figure 3 represents the overall computer program hierarchy. The main 
program establishes the link between subroutine (KREF), for the kinetic model 
of the reformer, and the following subroutines which determine the system 
performance and the mass and energy balance at various locations: BURN, CDPH, 

COMP, CON, CONV, DIVID, DMIX, ENFU, ENRE, ENSH, EQUK, FLAME, FUCE, HEPD, HEXC, 
PDFU, POSH, PUMP, PUP, REF, SNAE, SEPAR, and SHIFT. 

3.1 Main Program 

The main program performs the following functions: 

A. It reads the following input data: the thermophysical properties of 

methane, methanol, naphtha, water, oxygen, hydrogen, carbon monoxide, carbon 
dioxide, and nitrogen; data related to various components of the fuel cell 
power plant. 

B. For a given fuel, i.e., methane, methanol or naphtha, it carries out 

an itterative procedure to determine the thermodynamic state of the gas streams 
at various locations in the system, and to calculate the system efficiency and 
electric and heat energies output and other performance parameters. 

These calculations are carried out for two cases. In the first case, 
the kinetic effect on the reformer performance is considered to be negligible. 
In this case, the main program carries out these calculations without calling 
subroutine (ENRE). In the second case, the kinetic effect on the reformer 
performance is taken into consideration. For this case, the main program calls 
subroutine (KREF) and bypasses subroutines (ENRE), (EQUK), (REF), and (SNAE). 
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C. It creates a printout of the input data, the results of thermodynamic 
states of the gas streams, the system performance parameters, the output heat 
and electric energies. 

The nomenclature for the main program is shown in Table 1, and the flow 
chart appears in Figure 4. 

The equations contained in the main program are given below: 

1. Calculate inlet air flow rate in the burner: 

DNSS(33,2) = (1+EXT*0.01)*(DNSS(14,3)+DNSS(14, 5 ) ) /2+CK*DNSS( 14, 1 ) ) (3-1-1) 

where CK = stoichiometric number of oxygen used to burn the fuel: 
for methane, CK = 2 
for methanol, CK = 1.5 
for naphtha, CK = 15 

2. Calculate the saturation pressure of water for a given temperature 
T( 22) = -B/ (AL0G( (DNSS( 1, 1 )*SMRA-DNSS(21 ,6) ) / (DNSS ( 1 , 1)*SMRA- 

TKNSS(21) )*P0PS)-A) (3-1-2) 

where A and B are constants which have the following values for water: 

A = 13.954316, atm 

B = 5204.9597, atm-K 

3. Calculate the output AC power for a given DC power 

AC = (-1 . 0148+SQRT ( 1 .0148*2-4*0 - 056/ 108* ( 0 . 0472*108-WK) ) ) / 2*0 . 0456/ 108) (3-1-3) 


4. Calculate the flow rate of cooling water used in condenser 

DNSS(36,6) = QQT ( 5 ) / 1 / 18/ (355-TAT) (3-1-4) 



25 


A 

AA1 

AA2 

AA3 

AHLU 

AHRN 

AIRL 

APPD 

ATMP 

B 

BPNA 

BSPAC 

CD 

CLENH 

CLEPD 

CLH 

CN 

DG 

DHIN 

DHO 

DP 


TABLE 1 

MAIN PROGRAM NOMENCLATURE 


Constant for calculating saturated condition of water, atm 
Thermal conductivity coeff. of gas I, Btu/hr-ft-R 
Viscosity coeff of gas I, Ibm/ft-hr 

Specific heat capacity coeff. of gas I. Btu/R-lb-mole of the form: 
AA3(1)+AA3(2)*T+AA3(3)*T 2 +AA3(4)/T 2 
Mole fraction of available hydrogen 
Percent free gas space 
Length of air channel, ft 

Total surface area of packing Acc. to the basis and oper. temp., ft 2 
Outlet temperature of gases, K 

Constant for calculating saturated condition of water, atm-K 

Boiling point of naphtha, C 

Baffle space, ft 

Current density. A/cm 2 

Length of tube in heat exchanger, ft 

Length of shift converter (0K=1), reformer (OK-2 for methanol and 
naphtha), ft 

Clearance in heat exchanger, ft 

U*A/CMIN in heat exchanger 

Standard free energy change, Cal/g-mole 

Enthalphy change due to temperature change of inlet fluid, Cal/g-mole 
Integration constant to calculate H 
Catalyst pellet diameter, ft 
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DPD 

DSHO 

DSN 

DSO 

DTH 

DX1 

DX2 

DX3 

DZZ 

EA 

EPS 

ERR 

EXA 

EXT 

FCO 

FLOAR 

FULE 

HNA 

I 


TABLE 1 

MAIN PROGRAM NOMENCLATURE 
(cont'd) 

Diameter of shift converter (JK=1), reformer (JK=2 for methanol and 
naphtha), ft 

Cathode inlet water of fuel cell, g-mole/hr 

Cathode inlet nitrogen of fuel cell, g-mole/hr 

Cathode inlet oxygen of fuel cell, g-mole/hr 

Fraction of Delta T over inlet gas film in the heat exchanger 

Outside diameter of reformer center tube, ft 

Inside diameter of outside reformer tube, ft 

Outside diameter of outside reformer tube, ft 

Increment height of finite difference model in the reformer, ft 

Activation enegy for Arrhenius expression, Cal/g-mole CH4 

Reactor void fraction 

Convergence criteria 

Fraction of extra air in fuel cell 

Fraction of extra air in burner 

Mole fraction of co contain 

Flow area in heat exchanger, 

Length of fuel channel, ft 
Specific heat of naphtha, Btu/lbm-R 
Gas number 

I = 1 Fuel (methane, methanol, naphtha) 

I = 2 Oxygen 
I = 3 Carbon Monoxide 
I = 4 Carbon Dioxide 
I = 5 Hydrogen 
I = 6 Water 
I = 7 Nitrogen 



TABLE 1 


MAIN PROGRAM NOMENCLATURE 
(cont'd) 

ISSH : ID of shell in heat exchanger, ft 

IDTH : ID of tube in heat exchanger, ft 

IFUEL : Fuel Type 

1 = Methane CH4 

2 = Methanol CH30H 

3 = Naphtha C7H16 

IDNO : Number of tri al-and-error loops 

IHUI : Stoichiometric number 

IP : Index of operation condition in the reformer and shift converters 

IP = 1 Adiabatic Operation 
IP = 2 Isothermal Operation 

KO : Frequency factor for Arrhenius expression, lb-mole CH^ /lb 

cata.-hr-atm 

NN : Stream number of exit of shift converter 

NOR : Scale factor in the model of reformer 

NPFU : Number of cell plates in the fuel cell stacks 

NPH : Number of tube passes 

NRH : Number of rows for tubes 

NTAA : Number of air flow channel in one cell plate 

NTAF : Number of fuel flow channel in one cell plate 

NTPD : Number of tubes in shift converter (JK=1), Reformer (JK=2 for 
methanol and naphtha) 

ODTH : OD of tube, ft 

OU : 02 utilization 

PAT : Ambient pressure, atm 
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TABLE 1 


PIN 

PINFU 

PITCH 

POP 

POUT 

PL 

RHOB 

S 

SITS2 

SK 

SKI 

SMRA 

SURFC 

SV(1) 

svw 

TACOA 

TACOF 

TAT 

TC 

TCAS 

TCBS 


MAIN PROGRAM NOMENCLATURE 



Inlet pressure, atm 

Inlet pressure of fuel cell stacks, atm 
Pitch of heat exchanger, ft 
Operation pressure, atm 
Outlet pressure, atm 
Platinum catalyst loading, mgPT/CM 2 
Bulk density of cata., lbs/ft^ 

Side length of an assumed square flow duct for combustion gas, ft 

Ratio of total inside-tube cross-sectional area per pass to header 
cross-sectional area per pass 

Equilibrium constant 

Equilibrium constant with pressure different from 1 atm 

Steam/fuel ratio 

Surface per line, ft 

Specific volume of fuel 1, ft^/lbm 

Specific volume of water, ft^/lbm 

Inlet air temperature of fuel cell stack, K 

Inlet fuel temperature of fuel cell stack, K 

Ambient temperature, K 

Critical temperature, K 

Total heat capacity constant A 

Total heat capacity constant B 

Total heat capacity constant C 


TCCS 



TABLE 1 


TDNS : 

TIN : 

TOP : 

TOVO : 

TOUT : 

VHNA : 

WAT : 

WIDAA : 

WIDAF : 

X : 

ZH : 

DINSC(I) : 

DNS(I) : 

HA{ J) : 

HCAS(I), 
HCBS(I), 
HCCS(I) : 

HS( I ) : 

NNS(I) : 

WM( I ) : 

DNSS( I , J) : 


MAIN PROGRAM NOMENCLATURE 
(cont 'd) 

total amount of material, g-mole 

Inlet fluid temperature, K 

Operation temperature, K 

Total volume of inlet flow, m 2 

Outlet temperature, K 

Vaporized heat of naphtha, Cal/g-mole 

Relative humidity of air, g water/g air 

Width of square air channel in the fuel cell stack, ft 

Width of square fuel channel in the fuel cell stack, ft 

Necessary amount of oxygen in cathode, g-mole/hr 

Reformer length, ft 

Inlet amount of gas I, g-mole 

Inlet (outlet) amount of gas I, g-mole 

Surface area of heat exchanger J, m 2 


Heat capacity const, of gas I, Cal/g-mole-K of the form 
HCAS+HCBS*T+HCCS*T 2 

Heat of formation of gas I at 298 K, 1 atm, Cal/g-mole 
Stoichiometric coefficient of gas I 
Molecular weight of gas I, g/g-mole 
Flow rate of gas J in stream I, g-mole/hr 
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* stream number (refer to Figure 1) 

Figure 4 Flow chart of executive program for simulating CSU's FAFC system 

steady state performance 






























Figure 4 continued 
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3.2 Subroutines 

BURN, CDPH , COMP, COND, CONV, DIVID, DMIX, ENFU, ENRE, ENSH, EQUK, FLAME, 
FUCE, HEPD, HEXC, PDFU, PDSH, PUMP, PUP, REF, SNAE, SEPAR, and SHIFT. 

A. Subroutine BURN: This subroutine calculates the mass balance across 

the burner. It is assumed that combustion goes to completion and that the 
anode exhaust fuels the burner with 200 percent stoichiometric air. The 
illustrated equations containes in BURN for methane input fuel are: 

1. Calculate the amount of oxygen reacted: 

X = 0.5*DNS(3)+0.5*DNS(5)+2*DNS(1) (3-2-1) 

2. Calculate the amount of carbon dioxide produced 

XY = DNS(3)+DNS(i) (3-2-2) 

3. Calculate the amount of water produced 

Y + DNS(5)+2*DNS(1). (3-2-3) 

4. Calculate the exit composition 

DNS(l) = 0 
DNS(3) = 0 
DNS(5) = 0 
DNS(2) = DNS(2)-X 
DNS(4) = DNS(4)+XY 
DNS(6) = DNS(6)+Y 
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Figure 5 Flow Chart of BURK 
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B. Subroutine CDPH: This subroutine calculates the heat transfer rate in 

the evaporator E-10 and in the condenser E-7. The equations contained in CDPH 
are: 

1. Calculate heat transfer rate in heat exchangers E-7 and E-10 

QT = ( ( I-TC1) / (1-0.577) )**0. 38*9700 

*DNSC(6)+(TCB-TC1)*1*18*DNSC(6) . (3-2-4) 

2. Calculate the boiling temperature of water at a given pressure 

TCB = ' B/ (A-ALOG(P) ) 

where A and B are constants referred to in Equation (3-1-2). 

C. Subroutine COMP: This routine calculate the power requirement and 

shaft work for the fuel compressor. The equations contained in COMP are: 

1. Calculate the compressor shaft work assuming adiabatic conditions 

WS = GAG*1.987*TIN*1.8*((£^)**((GAG-1)/GAG)-1)/(GAG-1) (3-2-5) 

2. Calculate the compressor shaft work assuming isothermal conditions 

WS = 1 . 987*T I N*1 . 8*AL0G ( POUT /PIN) (3-2-6) 

3. Calculate the compressor power requirements 


POW = WS*TDNS/641400 


(3-2-7) 




Figure 6 Flow Chart of CDPH 
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Figure 7 Flow Chart of COKP 
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D. Subroutine COND: This subroutine calculates the heat transfer duty in 

the condenser. The hot side stream is a gas mixture that contains steam. 

COND contains the following equations: 

1. Calculate the condenser heat transfer duty (sensible heat only) 

QT = QT+DNSH( I )*(HCAS( I )*(THI-THO)+HCBS( I )*(THI**2-TH0**2) 

+HCCS( I )*(THI**3-THO**3) . (3-2-8) 


2. Calculate the condenser heat transfer capacity with the Watson 
correction for latent heat 

QT = QT+( ( 1— ( THO/ 647. 1) ) / ( 1-0. 577)**0. 38*9700*DNSH (6) (3-2-9) 

where Watson correction is given as, 

<h f ->, 1-T 2 °' 38 

wSrf = ( t^tt7 ) (3 - 2 - 10) 


where : molar heat of vaporization at condition i 


T ri ’ 


reduced temperature at condition i. 


E. Subroutine CONV: This subroutine finds the roots of the nonlinear 

equation x=f(x) by the Wegstein iteration scheme which accelerates convergence 
to the roots provided f(x) has a continuous first derivative. CONV contains 
the following equation: 

1. Calculate the roots of a given nonlinear function: 

XT = (XA( NR)*YV-YA(NR)*XV )/ (XA( NR)-XV + YV-YA(NR) ) (3-2-11) 


F. Subroutine D I V I D : This subroutine calculates the material balance 

around the divider with known divider factor. It is assumed that there is no 
temperature change in the streams and that specific enthalpy remains constant. 
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Figure 8 Flow Chart of COMD 
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Figure 9 Flow Chart of DKDC 
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G. Subroutine DMIX: This subroutine calculates mass and energy balances 

around the mixer through which two streams combine to produce a single 
stream. The flow chart for DMIX is shown in Figure 8. DMIX contains the 
following equations: 

1. Calculate the outlet temperature 

TOUTC = (TCASl*TINl+TCAS2*(TIN2-T0UT)+TCBSl/2.* 
(TIN1**2-T0UT**2)+TCBS2/2.*(TIN2**2-T0UT**2)+TCCS1* (3-2-12) 

(TINl**3-TOUT**3)/3.+TCCS2*(TIN2**3-TOUT**3)/3.)/TCASl 

2. Calculate the outlet pressure 

POUT = (TDNS1+TDNS2) / (TDNS1*TINI /PINI + TDNS2*TIN2/PIN2)*T0UT (3-2-13) 

H. Subroutine ENFU: This subroutine uses mass and energy balances in the 

fuel cell to calculate the following performance parameters: operating 

voltage, open circuit voltage, free energy change at fuel cell operating 
conditions, heat of reaction for methane, heat of reaction for methanol, heat 
of reaction for naphtha, fuel cell outlet temperature and stream composition, 
electrical work produced, heat energy rejected, voltage efficiency, thermo- 
dynamic efficiency, heating value efficiency, and fuel cell efficiency. The 
flow chart for ENFU appears in Figure 10. Subroutine ENFU contains the 
Equations (2-2-5) to (2-2-18). 

I. Subroutine ENRE: This subroutine is used to calculate the energy 

balance of reformer in lumped model. This model was based on the assumption 
that all chemical reactions reach equilibrium at the input temperature 
(isothermal operation) or the average temperature (adiabatic operation). Then 




Figure 10 Flow Chart of ENFU 
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Figure 10 continued 
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the energy balance contains the sensible enthalpy change and the enthalpy 
change of reactions. The mathematical model was described in the Equations 
(2-1-10) to (2-1-12). The flow chart of ENRE is shown in Figure 11. 

J. Subroutine ENSH: This subroutine is used to calculate the energy 

balance of shift converters (both high temperature converter and low 
temperature converter). Since methanol fuel does not need the shift 
converter, this subroutine will be skipped when input fuel is methanol. The 
mathematical model and flow chart of ENSH are shown in the Equation (2-1-6) 
and Figure 12, respectively. 

K. Subroutine EQUK: This subroutine calculates the equilibrium constants 

of the process gases in the demethanation. and water shift reactions. The 
mathematical model for the equilibrium constant was based on the Van't Hoff 
equation 

d 3n K = dT (3-2-14) 

Rr 

This equation can be integrated after expressing aH° in terms of the specific 
heats of the stream gases to yield, 

5n K = ^ Jln(T) + T + T 2 + AI (3-2-15) 

where ab and A-y are total heat capacity constants in the specific reaction, 

DHO and AI are constants of integration which can be evaluated from the 
standard enthalpy and standard free energy change. The flow chart for EQUK 
appears in Figure 13. The equations contained in EQUK are: 




Figure 11 Flow Chart of ENRE 
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Figure 12 Flow Chart of MSH 
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Figure 13 Flow Chart of EQUK 
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1. Calculate heat capacity constants 

TCAS = TCAS+NNS( I )*HCAS( I ) 

TCBS = T CBS+NNS ( I )*HCBS( I ) (3-2-16) 

TCCS = TCCS+NNS( I )*HCCS( I ) 

2. Calculate enthalpy of reaction change 

DH = DH+NNS( I )*HS( I ) (3-2-17) 

3. Calculate free energy of reaction 

DG = DG+NNS( I )*GS( 1 ) (3-2-18) 

4. Calculate constant DHO 

DHO = DH-TCAS*TST-TCBS*TST**2/2-TCCS*TST*3/3 (3-2-19) 

5. Calculate constant AI 

AI = ( DHO-DG-TCAS*TST*ALOG(TST )-TCBS/2*TST**2 

-TCCS/6*TST**3) /TST/R (3-2-20) 

6. Calculate equilibrium constant 

SK = EXP(-DH0/R/T0P+TCAS/R*AL0G(T0P)+TCBS/2*T0P/R 

+TCCS/6/R*T0P**2+AI) (3-2-21) 

L. Subroutine FLAME: This subroutine calculates the sensible enthalpy, 

the enthalpy change of reaction, and the maximum flame temperature in the 
burner. In the derivation of the mathematical model, it was assumed that the 
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combustion process goes to completion with negligible dissociation of the 
products and 200 percent stoichiometric air. The flow chart for FLAME appears 
in Figure 14. FLAME contains the following equations: 

1. Calculate the enthalpy of reaction change at 298 K 

DH = DH+DNS( I )*HS( I ) -DINS ( I )*HS( I ) (3-2-22) 

2. Calculate the sensible enthalpy change 

DH = DH+DINS( I )*(HCAS(I )* (298-TIN) + HCBS( I ) /2*( (298) 

**2-TIN**2)+HCCS(I )/3*( (298)**3-TIN**3) ) (3-2-23) 


3. Calculate total heat capacity constants 
TCAS = TCAS+DNS(I)*HCAS(I) 

TCBS = TCBS+DNS(I)*HCBS(I) (3-2-24) 

TCCS = TCCS+DNS ( I ) *HCCS ( I ) 


4. Calculate the adiabatic FLAME temperature 

TFC = (-DH-TCBS/2*( (TF)**2-(298)**2)-TCCS/3* 

( (TF )**3-(298)**3) ) /TCAS+298 (3-2-25) 


M. Subroutine FUCE: This subroutine calculates the mass balance in the 

fuel cell stack, which is described in Equations (2-2-2) to (2-2-4). Flow 
chart of FUCE is shown in Figure 15. 

N. Subroutine HEPD: This subroutine calculates the pressure drop in the 

heat exchangers used in the fuel cell power plants. It was assumed that BWG14 
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Figure 14 Flow Chart of FLAKE 





figure 15 Plow Chart of FUG 
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tubes with nominal size of 3/4 inch were used in the heat exchangers. The flow 
chart for HEPD appears in Figure 16. HEPD contains the following equations: 


1. Calculate the number of tubes 

NT = HA/Q. 3048**2 /NP/CLEN /SURF C (3-2-26) 

2. Calculate number of baffles 

NB = CLEN/BSPAC (3-2-27) 

3. Calculate free area between baffles 

FAREA = IDS/ (ODT+CL )*CL*BSPAC (3-2-28) 

4. Calculate ratio of pitch, transverse to flow, to tube diameter 

XT = PITCH/ODT (3-2-29) 

5. Calculate friction factor 

FPRI = SB0*(0DT*GS/AMUI )**(-0. 15) (3-2-30) 

6. Calculate pressure drop 

DP = B0*2*FPRT*NR*GS**2/32 .174/ 3600**2 / RH 0/ 2116.2 (3-2-31) 


0. Subroutine HEXC: This subroutine calculates the energy analysis in 

the parallel, counter and crossflow heat exchangers. From the assumption 
described in Section 2.1.1, the counter mode will be the only option used for 
heat exchangers in the system. Mathematical model was shown in the Equations 
(2-1-1) to (2-1-4). 
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Figure 16 Flow Chart of HEPD 
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P. Subroutine POFU: This subroutine calculates the pressure drops in the 

fuel channels and air channels. Dimensions of the fuel-cell stack are based 
on Westinghouse Stock No. 522. The mean pressure drop is evaluated by taking 
average of calculations based on inlet and outlet gas compositions. 

Q. Subroutine PDSH: This subroutine calculates the pressure drop in the 

packed reactors which are reformer and shift converters in our system. Ergun 
equation stated in Equation (2-1-7) was used to calculate pressure drop of 
reacting fluid caused by flowing through the packings. 

R. and S. Subroutines PUMP and PUP: Subroutine PUMP calculates the power 

required to pump water to a given pressure. PUP calculates the power required 
to pump naphtha or methanol to a given pressure. PUMP contains the following 
equation: 

1. Calculate the power required to pump water 

POW = SVW*144.*5.05051*0.0000001*WM(6)*14.7* 

(POUT— PIN)*DNS(6) /453.6 (3-2-32) 

T. and U. Subroutines REF and SNAE: Material balance in the reformer at 

the equilibrium state (lumped model) is analyzed in subroutine REF. Subroutine 
SNAE solves two nonlinear algebraic equations generated in REF. These two 
subroutines were more likely for the system with methanol or naphtha input 
fuel, whereas the kinetic model (Section 2. 1.3. 2) was used for the system with 
methane input fuel . 
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The material balances for methanol or naphtha input fuel in the reformer 
are similar to the discussion in the Section 2. 1.3.1, where methane input fuel 
was illustrated. 


Newton-Raphson method for solution of nonlinear algebraic equations is 
used in SNAE repeatedly to approach the equilibrium conversions of two 
parallel reactions (demethanation ana water shift reactions). The general 
description of Newton-Raphson method is as follows (Ref. 12), for two equa- 
tions f^Xp X 2 ) = 0 and f2 (Xp X g ) = 0: 


X 1 NEW = X 1 old + aX 1 
X 2 NEW = X 2 old + *X 2 


(3-2-33) 


where 


aX x = 


aX 2 = 


: !!] 
2 ax. 


- f. 


3x 2 


D 


afp 3f 1 

f - f i 

1 dx 1 2 ax 1 


(3-2-34) 


and D is the determinant of coefficient matrix (the Jocobian), which equals to 

d f i 3f 2 dfi df 2 

ax x 3x 2 ax 2 ax x 


V. Subroutine SEPAR: This subroutine calculates the outlet compositions 

in the liquid-vapor separator. The liquid-vapor equilibrium constant at given 
temperature is determined by Raoult's law which states 


XW = ( TDNS-DNS ( 6 ) / ( DK-1 ) 


(3-2-35) 
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( RETURN ) 


Figure 17: Flow Chart of Subroutine REF 







where XW: 
DK: 
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amount of water in liquid phase 

equilibrium constant of liquid-vapor system, which equals to 
(PSAT/POP) 


W. Subroutine SHIFT: SHIFT calculates the material balance in the 

shift converter. As discussed in Section 2.1.2, one reaction, water shift 
reaction, dominates the material change in the shift converter. The 
mathematical model described in the Equation (2-1-5) will be solved by 
Newton's method (Ref. 12) in this subroutine. 

X. Subroutine V I NEW : This subroutine calculates the characteristic of 

current density and operating voltage in the PAFC stack. The cell voltage can 
be expressed as an explicit function of reactions, products, and current 
density (Section 2.2.2), while the calculation of the current density involves 
a trial and error procedure. The mathematical model is shown in Equations 
(2-2-8) and (2-2-13). 

3.3 Subroutine KREF 

The mathematical model developed in Section 2. 1.3. 2 was used to develop 
a Fortran computer code, which consists of an executive program (KREF), three 
subroutines and eleven functions. Finite difference method will be applied 
to solve these simultaneous differential equations (Equations (2-1-14) to 
(2-1-17)), with the inlet conditions as the boundary conditions. The defini- 
tion of finite difference section is expressed in Figure 18 and the summary of 
the basic difference equations is shown in Table 2, where the nomenclature of 
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variables and functions in the program is listed in Table 3. Figure 19 shows 
the flow chart of this program. 

KREF is similar to REPENT developed by Westinghouse (80-9E6-PAMEC-RI ) . 

A more detailed discussin of the subroutines and functions is given in Ref. 3 
and will not be repeated here. 
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TABLE 2 

SUMMARY OF THE BASIC EQUATIONS 
USED IN THE KINETIC MODEL OF THE REFORMER 


1. Demethanation Reaction Kinetic Mass Balance 

Xl( i+1) = XI ( i ) = -n— c~ K o _EA/R t TA ( 1+1 ) +460 l XMCOMP(i.l) 

u o o 

2. Water Gas Shift Equilibrium 

_ -B- B 2 -4CA 
~ 2A 

with A = [K2( i+l)-l][F3+Xl( i+1) Fl]2 

B = [F3+Xl(i+1) F1][2X1( i + l) FI K2( i+1 )-F2 K2(i+1)-F4 K2( i+l)-5Xl( i+1) 
F1-F3-F5] 

C = K2( i+1) F2 F4-2K2 ( i+1 ) F2 Xl( i+1) F l-[ F3+X1 ( i +1 ) ] [F5+4X1( i+1) FI] 

3. Reforming Gas Energy Balance 

TC( i+1) = T - H - ^ +1 ^ [AM-MH CH] + [AM+MH CH] - TC(i) 

with AM = (UI tt D2 aZ)/2 

4. Combustion Gas Energy Balance 

TH(m , =TH(i) AN^m. TC(j) im.AL 

wherelFACP: sum of the component's heat capacity in the reforming gas 

AL = FI ( — DH1 )[Xl(i+l )— XI ( i )]+[F3+Xl(i+l) FI] (-DH2) [X2( i+1 )-X2( i ) ] 

AN = MH C iJl EFA - CP - ) - -LFACP+MH CH 


5. Pressure Drop of Reforming Gas 

P ( i+1 ) = P ( i )-aP 

where a P can be calculated from Ergun equation (Equation (2-1-7)). 
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TABLE 3 

NOMENCLATURE OF VARIABLES AND FUNCTIONS IN SUBROUTINE KREF 


Component 

2 

3 

4 

5 

6 
7 


Compound 

CH4 

CO 


C0 2 

h 2 o 


0 2 


Variable 


Definition 


FL for J = 0 to 7 Total reformer gas feed flow rate and component 

feed flow rates for components 1 thorugh 7 
(molar basis) 


MH 


Total combustion gas flow rate (molar basis) 


CGCOM for 0 = 1 to 7 

K1 

K2 

KO 

WM 

EA 

RHOB 

EPS 

D1 

D2 

D3 

S 

DP 

P 


Combustion gas component flow rates for compo- 
nents 1 to 7 (molar basis) 

Equilibrium constant for demethanation reaction 

Equilibrium constant for water shift reaction 

Frequency factor for Arrhenius expression 
k + k Q exp(-Eact/RT) 

Molecular weight 

Activation energy for demethanation reaction 

Catalyst bulk density 

Reactor void volume 

Reformer center tube outside diameter 

Reformer outer tube inside diameter 

Reformer outer tube outside diameter 

Characteristic dimension of the combustion gas 
flow duct (geometry is square) 

Catalyst particle oiameter 

Pressure 
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TABLE 3 

NOMENCLATURE OF VARIABLES AND FUNCTIONS IN SUBROUTINE KREF 



(cont 'd) 

Variable 

Definition 

TCO 

Reformer gas feed temperature 

THZ 

Combustion gas feed temperature 

Z 

Length of reformer tube 

X1(I) 

Kinetic conversion by demethanation reaction in 
Increment I 

XE2( I ) 

Actual conversion by water shift reaction in 
Increment I 

CO 

Initial methane concentration in reformer gas 

UO 

Initial reformer gas linear velocity 

T 

Temperature 

TK2 

Equilibrium constant for demethanation reaction 
at temperature T 

TX2 

Equilibrium conversion for water shift reaction 
at temperature T 

TDH1 

Heat of reaction for demethanation reaction at 
temperature T 

TDH2 

Heat of reaction for water shift reaction at 
temperature T 

X2 

Equilibrium conversion for water shift reaction 

DH1 

Heat of reaction for demethanation reaction 

DH2 

Heat of reaction for water shift reaction 

TF 

Total moles of reformer gas 

TVIS 

Viscosity at temperature T 

VIS 

Viscosity 

THC 

Thermal conductivity 
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TABLE 3 

NOMENCLATURE OF 

VARIABLES AND FUNCTIONS IN SUBROUTINE KREF 


(cont'd) 

Variable 

Definition 

TTHC 

Thermal conductivity at temperature T 

HI 

Inside heat transfer coefficient for the 
reformer outer tube 

THI 

Same as HI evaluated at temperature T 

HO 

Outside heat transfer coefficient for the 
reformer gas outer tube 

THO 

HO evaluated at temperature T 

UI 

Overall heat transfer coefficient for reformer 
outer tube 

TUI 

UI evaluated at temperature T 

FCP 

Reformer gas heat capacity 

TFCP 

Reformer gas heat capacity at temperature T 

THO 

Combustion gas outlet temperature 

DZZ 

Incremental length 

TA( I ) 

Average temperature in Increment I 

I 

Increment counter 

TAK1 

K1 evaluated at TA(I, K) 

TAK2 

K2 evaluated at TA(I, K) 

TAX2 

Conversion X2 for water shift reaction assuming 
total system equilibrium 

XF 

Total number of moles in the reformer gas 

XU I 

Overall heat transfer in Incrememt I 

XCOMP( I , J) 

Moles of component J in Increment I 

COM(J) 

Feed moles of component J 
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TABLE 3 


NOMENCLATURE OF VARIABLES AND FUNCTIONS IN SUBROUTINE KREF 



Vari able 
XMCOMP ( I , J) 
TC( I ) 

TH( I ) 

TP(I) 

XDH1 

XDH2 

XV1S 

XTHC 

XHI 

XCGVIS 

XCGTHC 

XHO 

RE 

Error 


Definition 

Mole fraction components J, Increment I 

Reformer gas temperature in Increment I 

Combustion gas temperature in Increment I 

Iteractive variable for TC ( I ) 

Value of DH1 in Increment I 

Value of DH2 in Increment I 

Value of VIS in Increment I 

Value of THC in Increment I 

Value of HI in Increment I 

Combustion gas viscosity in Increment I 

Combustion gas thermal conductivity in 
Increment I 

Value of HO in Increment I 
Reynolds number 

Convergence criterion on methane conversion 
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Combustion 

Gases 


Combustion 

Gases 


XI (1+1 ) 
XE2(l+l) 

TC( 1+1 ) 
P(X+1) 


X1(I) 

XE2(I) 

TC(l) 

P(I) 



♦ 


th(i+i ) 


xui(l+1 ) 

th(i) 


Figure 18 Single Tube Kinetic Reformer Model 
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3.4 Program Operation 

The program input only consists of a NAMELIST data deck which must be in a 
specified order. The first NAMELIST set is called OPFC and contains the three 
values of the operating conditions in the fuel cell stack. The order of input 
data inside one NAMELIST need not be fixed. 

The second set ( I N IT ) contains the 11 values of the amount of the input 
fuel and ambient temperature and pressure. The dimension orders in the 
variables for the properties of gas mixture are fixed, which are (1) methane, 
(2) oxygen, (3) carbon monoxide, (4) carbon dioxide, (5) hydrogen, (6) water, 
and (7) nitrogen. 

The third set (CONDT) carries’ the information for the system operation, 
which includes the kind of input fuel, trial and error criterion, relative 
humidity, and extra percentage of needed air in the burner and the fuel cell 
stack. 

The fourth set (REPEN) contains the information for the kinetic model of 
reformer. These are the dimensions of reformer and the catalyst kinetic data 
used in the reformer. 

The fifth set (HEATX) contains the operating conditions for all the heat 
exchangers in the system and the transfer areas designed in the condenser and 


cooler. 
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The sixth, seventh, and eighth sets (HEPDC, PDSHH, and PDFUH) contain the 
dimensions of the heat exchanger, the shift converters, and the fuel cell 
stack, respectively. These data will be used to calculate the pressure drops 
in these three components. 

The last NAMELIST set (CATAI) specifies the kinetic data of the catalyst 
used in the fuel cell stack. 



INPUT DATA FOR SIMULATION OF CSU PAFC SYSTEM STEADY STATE PERFORMANCE 
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CONDT IFUEL 1 =1 Input fuel is methane 

= 2 Input fuel is methanol 
= 3 Input fuel is naphtha 

CONDT ERR O.Ol Criterion of convergence in system trial 

and error procedure 



TABLE 4 (cont'd) 

INPUT DATA FOR SIMULATION OF CSU PAFC SYSTEM STEADY STATE PERFORMANCE 

(SAMPLE RUN) 


NAMELIST 

LIST 

VARIABLE 

NAME 

SAMPLE 

DIMENSION VALUE 

UNIT 

DEFINITION 

CONDT 

IP 

2 


= 1 adiabatic operation in shift converters 
= 2 isothermal operation in shift converters 

CONDT 

I 

7 


Number of components in whole system 

CONDT 

EXT 

100 


Extra percentage of needed air in burner 

CONDT 

WAT 

0.015 

g water/g air 

Relative humidity of air 

CONDT 

EXA 

100 


Extra percentage of air in fuel cell stack 

REPEN 

ZH 

6 

ft 

Height of reformer 

REPEN 

DX1 

0 

ft 

Outside diameter of regenerative tube 

REPEN 

DX2 

0.15 

ft 

Inside diameter of reforming tube 

REPEN 

DX3 

0.1667 

ft 

Outside diameter of reforming tube 

REPEN 

KO 

10400 

lb-moleCH4/ 
lb cata-hr-atm 

Rate constant of demethanation reaction 

REPEN 

EA 

20000 

Cal/g-mole CH^ 

Activity energy of demethanation reaction 

REPEN 

RHOB 

80 

lb/ft 3 

Density of packing in refomer 

REPEN 

EPS 

0.487 


Void fraction in reformer 

REPEN 

S 

0.25 

ft 

Width of combustion gas square duct 

REPEN 

DP 

0.00328 

ft 

Diameter of catalyst in reformer 

REPEN 

DZZ 

0.25 

ft 

Height of finite-difference section 

HEATX 

CN 

1.3 

m 2 -K 

QxA/C min in heat exchanger 


cn 

to 



TABLE 4 (cont'd) 

INPUT DATA FOR SIMULATION OF CSU PAFC SYSTEM STEADY STATE PERFORMANCE 

(SAMPLE RUN) 


NAMELIST 

LIST 

VARIABLE 

NAME 

DIMENSION 

SAMPLE 

VALUE 

UNIT 

DEFINITION 

HEATX 

U 


48825.1 

cal / m^-hr-K 

o 

Overall heat transfer coefficient in heat 
exchanger 

HEATX 

HA 

7 

0.2 

L 

m 

Transfer area in E-7 



10 

0.2 

2 

m 

Transfer area in E-10 

HEPDC 

NPH 


2 


Number of tube passes 

HEPDC 

NRH 


5 


Number of tube rows 

HEPDC 

BASPC 


1 

ft 

Buffle space 

HEPDC 

ODTH 


0.0625 

ft 

O.D. of tube 

HEPDC 

PITCH 


0.0833 

ft 

Pitch of heat exchanger 

HEPDC 

CLH 


0.0208 

ft 

Clearance in heat exchanger 

HEPDC 

IDSH 


0.833 

ft 

I.D. of shell 

HEPDC 

IDTH 


0.04667 

ft 

I.D. of tube 

HEPDC 

FLOAR 


0.001716 

ft 

Flow area in heat exchanger 

HEPDC 

SURFC 


0.1466 

ft 

Surface area per line 

HEPDC 

CLENH 


2 

ft 

Length of tube 

HEPDC 

SITSZ 


0.5 


Ratio of total inside tube cross-sectional 
area per pass to header cross-sectional area 
per pass 

HEPDC 

DTH 


0.7 


Fraction of aT over inlet gas film in heat 
exchanger 



TABLE 4 (cont'd) 

INPUT DATA FOR SIMULATION OF CSU PAFC SYSTEM STEADY STATE PERFORMANCE 

(SAMPLE RUN) 


NAMELIST 

VARIABLE 


SAMPLE 



LIST 

NAME 

DIMENSION 

VALUE 

UNIT 

DEFINITION 

PDSHH 

DPD 

1 

1.18 

ft 

Diameter of shift converters 

PDSHH 

AHRN 

1 

0.66 


Void fraction in shift converters 

PDSHH 

APPD 

1 

69 

ft 2 

Total surface area of packing 

PDSHH 

CLEPD 

1 

5.91 

ft 

Length of shift converters 

PDSHH 

NTPD 

1 

1 


Number of tubes in shift converters 

PDFUH 

NTAF 


140 


Number of fuel flow channels in stack 

PDFUH 

FULE 


1.42 

ft 

Length of fuel channel 

PDFUH 

WIDAF 


0.00974 

ft 

Width of square. fuel channel 

PDFUH 

NPFU 


3365 


Number of cell plates 

PDFUH 

NTAA 


40 


Number of process air flow channels 

PDFUH 

AIRL 


1 

ft 

Length of air channel 

PDFUH 

WIDAA 


0.00515 

ft 

Width of square process air channel 

CATAI 

SRO 


0.44 

^-cm 

Cell resistance at 450° K 

CATAI 

SA 


400 

cm 2 /mg 

Surface area of catalyst 

CATAI 

CU 


0.15 


Utilization of catalyst 

CATAI 

CL 


0.75 

mg /cm 2 

Catalyst loading 

CATAI 

ALFA 


0.50 


Transfer coefficient 

CATAI 

SN 


2 

g-equi valent 

Number of Faraday equivalents transferred 

CATAI 

FCONST 


96500 

C/g-equi valent 

Faraday constant 

CATAI 

DKC 


240000 

A/ atm 

Constant to calculate limiting current 
density 
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All of the input variables are listed in Table 4, along with their units 
and numerical values in the sample run, which will be discussed in the 
following section. 

3.5 Sample Problem 

The computer codes developed in previous sections were combined to 
simulate the PAFC system performance. The lumped model of each component was 
used to simulate the CSU design (Figure 1), except in the reformer where the 
distributed model was used for the methane input fuel. 

This sample problem was to simulate the 110 kW AC PAFC system with methane 
input fuel. The input data, which is discussed in the previous section, is 
displayed in Figure 20. Figure 21 contains the output generated by the sample 
data input. First, the program reprints all of the input data. Next, the 
program prints out the operating conditions (temperature and pressure) of 
reformer, shift converters, and liquid separator. For the fuel cell stack, 
the printout will contain the operating temperature, operating pressure, open 
circuit potential, operating potential, current density, catalyst loading, 
fuel and oxident utilizations, the different kinds of efficiency, DC and AC 
electrical work, and heat released from the stack. 

Next, the P-T-V (V as a flow rate) status of each stream numbered in 
Figure 1 will be listed on a new printout page. The last piece of information 
printed is the duty, transfer area, and efficiency of each heat exchanger 
numbered on flow diagram, and power spent in the air and fuel compressors and 


pump. 
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SOPFC T0PFC=443. ,UT=0.8,CD=325. , 

SEND 

SIN IT DNSM=1216 .,0.,1.360,21.8,166.,0.,0. , TAT = 298. ,PAT=1 . , SMRA = 3 . , P0PR=5 . 0 
SEND 

SCONDT I FUEL = 1, ERR = 0 . 0 1 , IP = 2 , I =7 , EXT = 1 0 0 . , WAT = 0 . 0 1 5 , EXA = I 0 0 . , 

SEND 

SREPEN ZH = 6 . ,DX1=0 . ,DX2=0 . 15,DX3 = 0 . 1667 , KO = 1 . 040E+04,EA=20000 . ,RHOB=80 . 

, EPS : 0. 487, S : 0. 25, DP : 0. 00328, DZZ = 0 . 25 , 

SEND 

SHEATX CN=1.3,U=48825.1,HA(7)=0.2,HA(10)=0.2, 

SEND 

SHEPDC NPH=2,NRH=5,BSPAC=1. , ODTH= . 0 6 25 , P I TCH= . 0833 , CLH= . 0208 , IDSH= . 833 , 
IDTH=. 04667, FLOAR=. 0017 16, SURFC=. 1466, CL ENH=2. , S 1TS2=0 . 5 , DTH=0 . 7 , 

SEND 

SPDSHH DPD=1 . 18 , 0 . , AHRN = 0.66, . 0 , APPD=6 9 . , 0 . , CLEPD=5 . 91 , 0 . , 

NTPD=1, 0 , 

SEND 

SPDFUH NTAF=140, FULE=1 .42,WIDAF=. 0 0 9744 , NPFU=3 36 5 , 

NTAA=40,AIRL=1. ,WIDAA=. 00515 
SEND 

SCAT AI SR0=.44, SA=400 . , CU= . 1 5 , CL = . 7 5 , AL FA = . 5 , SN = 2 . , FCONST = 96500 . , 

DKC=2 . 4E5 , 

SEND 


Figure 20 Sair.pl e Input Data 
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THE STEAM/METHANE RATIO IN THE REFORMER IS 3.00 

THE REFORMER IS OPERATING UNDER THESE CONDITIONS 
INLET PRESSURE 9.90ATM OUTLET PRESSURE 9.72ATM 
INLET TEMP. :1035.87 K OUTLET TEMP. : 920.63 K 

THE HIGH TEMP. SHIFT CONVERTER IS OPERATING UNDER THESE CONDITIONS 

OPERATING TEMP.: 755.26 K 

OPERATING PRESSURE: 9.72ATM 

INLET TEMP . : 755.26 K 

OUTLET TEMP. : 755.26 K 

THE LOU TEMP. SHIFT CONVERTER IS OPERATING UNDER THESE CONDITIONS 

OPERATING TEMP.: 399.68 K 

OPERATING PRESSURE: 9.72ATM 

INLET TEMP. : 399.68 K 

OUTLET TEMP. : 399.68 K 

THE LIQUID SEPERATER IS OPERATING UNDER THESE CONDITIONS 
OPERATING TEMP.: 357.39 K 
OPERATING PRESSURE: 9.92ATM 

THE FRACTION OF CO IN THE FEED IS 0.00029 

THE FUEL CELL IS OPERATING UNDER THESE CONDITIONS 

THE OPERATING TEMPERATURE : 993. 00K 

THE OPERATING PRESSURE: 4.79ATM 

THE OPEN CIRCUIT POTENTIAL: 1.171 V 

THE OPERATING POTENTIAL: 0.658 V 

THE CURRENT DENSITY: 0 . 325 A/CMX*2 

THE CATALYST LOADING: 0 . 750PT/CM*«2 

THE FUEL UTILIZATIONS. 800 

THE OXYGEN UT I L IZAT ION : 0 . 500 

THE ANODE SIDE INLET TEMP. IS 399.68 K 

THE CATHODE SIDE INLET TEMP. IS 386.91 K 

THE THERMODYNAMIC EFFICIENCY OF FUEL CELL IS 0.92897E OOTHE CURRENT EFFICIENCY IS 0.80000E 00 

THE VOLTAGE EFFICIENCY IS 0.56251E OOTHE HEATING VALUE EFFICIENCY IS 0.90970E 00 

THE FUEL CELL EFFICIENCY ISO. 3780 

THE ELECTRICAL WORK IS 0.12601E 03KW 

THE TOTAL HEAT RELEASE IS 0.53266E 08CAL 


THE AC OUTPUT IS 113.77KU 


Figure 21 continued 
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THE DUTY OF HEAT EXCHAHGER (CAL.) (0 MEANS NO THIS NO. HEAT EXCHANGER) 

E- 1 E-2 E-3 E-A E-5 E-6 E-7 E-8 E-9 E-10 

0.A5316E 07 0.99787E 07 0.1A123E 08 0.33156E 08 0.10910E 09 0.10777E 08 0.53266E 08 0.00000 0.51060E 07 0.37070E 08 

THE SURFACE AREA OF HEAT EXCHANGER (M**2) (0 MEANS NO THIS NO. HEAT EXCHANGER) 

E- 1 E-2 E-3 E-A E-5 E-6 E-7 E-8 E-9 E-10 

0.38308E 00 0.135A1E 01 0.1A96AE 01 0.22603E 01 0.21182E 02 0.1A996E 01 0.20000E 00 0.00000 0.82A20E 00 0.20000E 00 

THE EFFICIENCY OF HEAT EXCHANGER (0 MEANS NO THIS NO. HEAT EXCHANGER OR IS CONDENSER) 

E- 1 E-2 E-3 E-A E-5 E-6 E-7 E-8 E-9 E-10 

0.68881E 00 0 . 59n*»6E 00 0.63181E 00 0.57A12E 00 0.00000 0.65302E 00 0.00000 0.00000 0.6A153E 00 0.00000 

THE POWER OF AIR COMPRESSOR: 6A . OAHP 

THE POWER OF METHANE COMPRESSOR: 3.80HP 

THE POWER OF PUMP : 0.001A6HP 


Figure 21 continued 
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3.6 Discussion 

There has been much interest in the effect of alternate commercial fuels 
on the performance and costs of PAFC power plant. The computer program 
developed can allow for methanol or naphtha as input fuel. The system with 
methanol input fuel obtains the highest efficiency among the three fuels, 
where 40-45 percent of the PAFC stack compared to 35-40 percent of methane 
input fuel. More detailed discussion on this topic is presented in Ref. 14. 

Since there are a lot of tri al-and-error procedures in the program, the 
infeasible initial guesses will cause the calculations looping or overflowing 

For the naphtha input fuel, because of computation problem (overflow) in 
the subroutine SNAE, the conversion will be assumed in the reformer. The 
problem will be amended as required. 

3.7 Further Developments 

This PAFC system steady state simulation program can be modified to allow 
different flow diagram. For example, it has been used to simulate the 
Westinghouse 7.5 MW PAFC power plant (Ref. 13), and the results were shown in 
Ref. 11. 


Further developments have been completed, which include 3-D temperature 
and current density distributions (distributed model) of fuel cel stack (Refs 
11, 15), kinetic model for regenerated-type reformer, distributed simulation 
of PAFC system steady state performance (Ref. 11), and the simulation of PAFC 
power plant system transient responses in the load changing period (Ref. 11). 
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Listing of the Steady State Performance Lumped Model 


0000020C BLOCK DATA 

00C0C40CC 

0000C60 C******BLOCK DATA FOR CSU PROGM FUEL CELL POKER PLANT PERFORMANCE 
OOCOIOOC DIMENSION GS< 7 ) ,HS( 7 1 ,HCAS(7) ,HC6S( 7 ) ,HCCS( 7 ) , WMt 7 ) , SV( 3 ) , 

0000120C 1HLHV(7),A1(2,7),A2(2,7),A3<4,7) 

000C140C COMMON /ETHDA/ G3,HS,HCAS,HC6S,HCCS 

0000160C COMMON /TO/ TC/CCNS/ A,B 

0000130C COMMON /GAG/ GAGM.GAGA 

000020CC COMMON /NAFH/HNA , BPNA > KMNA i VHNA 

0000220C COMMON /UM/ KM 

0000240C CCMMCN/SV/ SV.SVW 

0C00260C COMMON/H LH VI/ HLHV 

0000230C COMMON/THCC/A1/VIFC/A2/HTCPC/A3 

000030CCC 

0000320CC STANDARD VALUE OF FREE ENERGY, ENTHALPHY AND HEAT CAPACITY CONSTANT 
0000340C DATA GS/-12140 ., 0 ., -32781 ., -94253 ., 0 ., -54635 0 . / 

0000360C DATA HS/-17839 0 -26 416 -94051 0 -57798 0 . / 

0C003S0C DATA HCAS/3 , 381 ,6 . 148 , 6 . 92 , 6 . 219 , 6 . 947, 7. 256 , 6 . 524/ 

0000400 DATA HCES/ . 018044 , . 0 03102 , . 001665 010396 0C02 , . 002298 00125/ 

0CCC420 DATA HCCS/-4 . 3E-06 , -9 . 2 3E -07 , -1 . 96E-07 , -3 . 545E-06 , A . 81E-07 , 2 . 83E-07 , 1 . E-09/ 

00G0440CC THE MOL. WEIGHT OF GASES 

Q000460C DATA KM/ 16 . ,32. ,28. ,40. ,2. ,18. ,28./ 

0000430CC THE THERMAL DATA OF NAPHTHA 

0C0050CC DATA HNA/0 .58/,EFNA/103. 4/, VKNA/7680 ./ 

C00052CCC THE CONST. FOR CAL. SATURATED PRESSURE <EXP(A-B/T)) 

C000540C DATA A , B/13 . 954316 , 5204 . 9597/ 

0000560CC CRITICAL TEMPERATURE OF WATER 

00005S0C DATA TC/646.447/ 

0C00600CC INSERT RATIO OF HEAT CAPASITY 
0000620C DATA GAGM/1 . 3/ , GAGA/1 . 4/ 

OOOC64CCC INSERT THE SPECIFIC VOLUME OF FUEL AND WATER 
0000660 DATA SV/0 .,. 02034, . 021/ ,SVW/ . 0162/ 

0000660CC INSERT THE HEATING VALUE 

00G07C0C DATA HLHV/- 191762 ., 0 ., -67636 ., 0 ., -57798 . ,0 ., 0 . / 

0000720CC INSERT THE SPECIFIC HEAT CAPACITIES, BTU/R-LB-MOLE 

OC0C740C OATA A3/5 . 34 ,6 . 39E-03 , 0 . , 0 . , 6 . 60 , 6 . 67E-C4 , 0 . , 0 . , 10 . 34 , 1 . 52E-03 , 0 . , 

0000760C 1-6 . 3342CE+C5,8 . 22 , 8 . 3E-05 ,4 . 136E-07 , 0 . ,6 . 62 ,4 . 5E-04 , 0 . , 0 . , 6 . 5 , 

0000780C 25.56E-04,0. ,0. , 6 . 732 , 8 . 36 E-03 , 5 . 53E -0 9 , 0 . / 

00008C0CC INSERT VISCOSITY, LEM/FT-HR 

0000820C DATA A2/3 . 4373E-05 , 2 . 5861E-02 , 3 . 8135E- 05 , 4 . 6688E-02 , 5 . 0368E-05 , 

0000840C 13.3634E-02,5.221E-05,1. 7532E-02 ,1 .8789E-05 ,2. 312E-02 ,4. 231SE-05, 

0000860C 24.6721E-02,2.45S0E-05,5.5613E-02/ 

C000880CC INSERT THERMAL CONDUCTIVITY , BTU/HR-FT-F 

00009C0C DATA Al/6 . 6539E-05 , 7 . 4614E-03 , 1 . 5349E-05 ,1 . 6581E-02 , 2 . 0451 E-05 , 

0000920C 11. 1726E-02, 4. 046 E-05, 1.8785 E-03, 1. 189 9E-04,1.0126E-01,1.6774E-05, 

0000940C 21 . 5182E-02 , 1.792E-05,1 -67E-02/ 

0000960C END 

,0000930 C 

0001000 C «*##********#»*»#»#**#***»#*#»#####**»*»#*»###***#***###******#**# 

0001020 C THIS PROGRAM IS USING FORTRAN TO SIMULATE THE PHOSPHORIC ACID FUEL 

0001C40 C CELL SYSTEM 

C 00 1060 C a#**#***##*#*#****#####*****###*######*****##*#******#***##******* 

0001080 C 

0001100 C DEFINITION: 
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0001120 C 
0001140 c 
0001160 C 
0001180 C 
0001200 C 
0001220 C 
0001240 C 
0001260 C 
0C01280 C 
0001300 C 
0001320 C 
0001340 C 
0001360 C 
0001330 C 
00014C0 C 
C001420 C 
0001440 C 
0001460 C 
CC01480 C 
0001500 C 
0001520 C 
0001540 C 
0001560 C 
00015S0 C 
0001600 C 
0001620 C 
0001640 C 
0001660 C 
0001630 C 
0001700 C 
0001720 C 
0001740 C 
0001760 C 
0001780 C 
0001800 C 
0001320 C 
0001840 C 
0001360 C 
0001830 C 
0001900 C 
0C01920 C 
0001940 C 
C001960 C 
0001930 C 
0C02000 C 
0002020 C 
0002040 C 
0002060 C 
C002CS0 C 
0002100 C 
C002120 C 
0002140 C 
C002160 C 
0002180 C 


A : CONSTANT FOR CAL. SATURATED CONDITION OF WATER 
AAl: THERMAL CONDUCTIVITY COEFF. OF GAS I, BTU/HR-FT-F 
AA2 : VISCOSITY COEFF. OF GAS I, LBM/FT -HR 
AA3: SPECIFIC HEAT CAPACITY COEFF. OF GAS I, BTU/R-LB-MOLE 
OF THE FORM : AA3! 1 )+AA3( 2 )#T+AA3( 3 )*T#*2+AA3( 4 )/T**2 
AHLU: MOLE FRACTION OF AVAILABLE HYDROGEN 
AHRN: PERCENT FREE-GAS SFACE 
AIRL: LENGTH OF AIR CHANNEL. FT 
ALFA: TRANSFER COEFF. 

APFD : TOTAL SURFACE AREA OF PACKING ACC. TO THE BASIS AND OPER. 

TEMPERATURE, FT**2 

ATMP: OUTLET TEMP. OF GASES, K 

B : CONSTANT FOR CAL. SATURATED CONDITION OF WATER 

BPNA: BOILING FOINT OF NAPHTHA, C 

BSPAC: BAFFLE SPACE, FT 

CD= CURRENT DENSITY ,MAMP/CM**2 

CL: CATALYST LOADING , MG/CM*#2 

CLENH: LENGTH OF TUBE IN HEAT EXCHANGER, FT 

CLEFD : LENGTH OF SHIFT CONVERTER! JK=1 ) ,REFORMER(JK=2 FOR METHANOL 
NAFHTHA ) , FT 

CLH : CLEARANCE IN HEAT EXCHANGER, FT 
CH: Q*A/CMIN IN HEAT EXCHANGER 
CU= CATALYST UTILIZATION 

DG-‘ STANDARD FREE ENERGY CHANGE .CAL/G-MOLE 

DH: STANDARD ENTHALFHY CHANGE AT REACTION , CAL/G-MOLE 

DHIN : ENTHALPHY CHANGE DUE TO TEMPERATURE CHANGE OF INLET FLUID 

CAL./G-MOLE 

DKC : CONST. TO CALC. LIMITING CURRENT DENSITY 
DHO: INTEGRATION CONSTANT IN CALCULATE H 
DP: CATALYST PELLET DIAMETER , FT 

DFD : DIA. OF SHIFT CONVERTER ( JK=1 ), REFORMER! JK=2 FOR METHANOL AND 
NAFHTHA), FT 

DSHO : CATHODE INLET WATER OF FUEL CELL, G-MOLE/HR. 

DSN: CATHODE INLET NITROGEN OF FUEL CELL, G-MOLE/HR. 

DSO: CATHODE INLET OXYGEN OF FUEL CELL, G-MOLE/HR. 

DTH : FRACTION OF DELTA T OVER INLET GAS FILM IN THE HEAT EXCHANGER 

DX1 : OUTSIDE DIAMETER OF REFORMER CENTER TUBE, FT 

DX2 : INSIDE DIAMETER OF OUTSIDE REFORMER TUBE, FT 

DX3: OUTSIDE DIAMETER OF OUTSIDE REFORMER TUEE, FT 

DZZ : INCREMENT HEIGHT OF FINITE DIFFERENCE MODEL IN THE REFORMER, 

EA: ACTIVATION ENERGY FOR ARRHENIUS EXPRESSION, CAL/GMOLE CH4 

EPS: REACTOR VOID FRACTION 

ERR - - CONVERGE CRITERIA 

EXA: FRACTION OF EXTRA AIR IN FUEL CELL 

EXT: FRACTION OF EXTRA AIR IN BURNER 

FCO : MOLE FRACTION OF CO CONTAIN 

FCCNST : FARADAY CONSTANT , 23061 CAL/VOLT-GM EQUIV. 

FLOAR : FLOW AREA IN HEAT EXCHANGER, FT**2 
FULE: LENGTH OF FUEL CHANNEL, FT 
HNA: SPECIFIC HEAT OF NAFHTHA, BTU/LBM-R 
I : GAS NUMBER 

1=1 FUEL) METHANE .METHANOL .NAPHTHA ) 

1=2 OXYGEN OR CARBON MONOXIDE! IN SUB. KREF AND 

RELATED SUBROUTINES) 
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0002200 C 1=3 CARBON MONOXIDE OR CARBON DIOXIDE! " ) 

0002220 C 1=4 CARBON DIOXIDE OR WATER! " ) 

00022*40 C 1 = 5 HYOROGEN 

00C2260 C 1=6 HATER OR NITROGEN! “ ) 

0002290 C 1 = 7 NITROGEN OR OXYGEN! '* ) 

0002300 C IOSH: ID OF SHELL IN HEAT EXCHANGER , FT 

0002320 C IDTH: ID OF TUBE IN HEAT EXCHANGER, FT 

0002340 C IFUEL- FUEL TYPE 

CC02360 C l: METHANE CH4 

C0C2380 C 2: METHANOL CH30H 

0032400 C 3= NAPHTHA C7H16 

0002*420 C IDNO • NO. OF REDO 

00024*40 C IHUi: STOICHIOMETRIC NUMBER 

0002460 C IP'- INDEX OF OPERATION CONDITION IN THE REFORMERINOT FOR METHANE) 
0002490 C SHIFT CONVERTOR 

0CC2500 C IP=1 ADIABATIC OPERATION 

0002520 C IP=2 ISOTHERMAL OPERATION 

0002540 C K0 = FREQUENCY FACTOR FOR ARRHENIUS EXPRESSION, LB MOLE CH4/LB 
0002560 C CATA.-HR-ATM 

0002530 C NNS STREAM NU-BER OF EXIT OF SHIFT CONVERTER 

0002600 C NCR: SCALE UP FACTOR IN THE MODEL OF REFORMER 

0002620 C NPFU: NO. OF CELL PLATE IN THE FUEL CELL STACK 

0002640 C NPH= NO. OF TUBE PASSES 

0002660 C NRH : NO. OF RCiiS FOR TUBES 

0002680 C NTAA: NO. OF AIR FLOW CHANNEL IN ONE CELL PLATE 

0002700 C NTAF: NO. OF FUEL FLOW CHANNEL IN ONE CELL PLATE 

0002720 C NT PD : NO. OF TUBES IN SHIFT CONVERTER! JK=1 ) .REFORMER! JK=2 FOR 

0002740 C METHANOL AND NAPHTHA) 

0002760 C ODTH: OD OF TUBE, FT 

0C02730 C CU: 02 UTILIZATION 

00C2800 C pat: AMBIENT PRESSURE, 1ATM 

0002620 C PIN: INLET PRESSURE, ATM 

0002640 C PINFU: INLET FRESSURE OF FUEL CELL STACK, ATM 

0002660 C PITCH: PITCH OF HEAT EXCHANGER, FT 

0002630 C FO?: OPERATION FRESSURE, ATM 

0002900 C FOUT: OUTLET FRESSURE, ATM 

0CC2920 C RHOB: BULK DENSITY OF CATA. , LB5/FT**3 

0002940 C S'- SIDE LENGTH OF AN ASSUMED SQUAPE FLOW DUCT FOR COMBUSTION GAS, 

0002960 C SITS2 : RATIO OF TOTAL INSIDE-TUBE CROSS-SECTIONAL AREA PER PASS TO 

0002930 C HEADER CROSS-SECTIONAL AREA PER PASS 

0003000 C SA: CATALYST SURFACE ,CM**2/MS 

0003020 C SK: EQUALIBRIUM CONSTANT 

0003040 C SKI: EQUALIBRIUM CONSTANT WITH PRESSURE DIFFERENT FROM 1 ATM 
C003060 C SMRA : STEAM/FUEL 

CC030S0 C SN: NUMBER OF FARADAY EQUIVALENTS TRANSFERED 
0003100 C SRO: CELL RESISTANCE AT 4S0 K,OHM-CM*»2 
0003120 C SUPFC: SURFACE PER LINE, FT 

OC 03140 C SV( I ) : SPECIFIC VOLUME OF FUEL I, FT**3/LBM 

‘0003160 C SVW : SPECIFIC VOLUME OF WATER, FT**3/LBM 

0003180 C TACOA: INLET AIR TEMP. OF FUEL CELL STACK, K 

C00320Q C TACOF: INLET FUEL TEMP. OF FUEL CELL STACK, K 

0003220 C TAT: AMBIENT TEMPERATURE, 293 K 

C003240 C TC= CRITICAL TEMPERATURE, K 

0C03260 C TCAS '• TOTAL HEAT CAPACITY CONSTANT A 
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0003230 C 

0003300 C 

0003320 C 

0003340 C 

0003360 C 

0003330 C 

0003400 C 

0003420 C 

0003440 C 

0003460 C 

0003480 C 

0003500 C 

C003520 C 

0003540 C 

0003560 C 

0003580 C 

0003600 C 

0C03620 C 

0C03640 C 

0003660 C 

0003680 C 

0003700 C 

0003720 C 

0003740 

0003760 

0003730 

0003300 

0003320 

OC 03840 

0003360 

0003330 

0003900 

0003920 

0003940 

0003960 

0003980 

0004000 

0004C20 

0004040 

0004060 

0004C80 

0004100 

0004120 

C004140 

C004160 C 

0004130 

00C4200 

0004220 

0004240 

0004260 

0004280 

0004300 

0004320 

0004340 


TCBS: TOTAL HEAT CAPACITY CONSTANT B 

TCCS: TOTAL HEAT CAPACITY CONSTANT C 

TONS: TOTAL AMOUNT OF MATERIAL. G-MOLE 

TIN: INLET FLUID TEMPERATURE, K 

TOP: OPERATION TEMPERATURE, K 

TOVO: TOTAL VOLUME OF INLET FLOW, M**3 

TOUT: OUTLET TEMPERATURE, K 

VKNA: VAPORIZED HEAT OF NAFHTHA, CAL/G-MOLE 

WAT : RELATIVE HUMIDITY OF AIR, G WATER/G AIR 

WIDAA : WIDTH OF SQUARE AIR CHANNEL IN THE FUEL CELL STACK, FT 

WIDAF : WIDTH OF SQUARE FUEL CHANNEL IN THE FUEL CELL STACK, FT 

X : THE NECESSARY AMOUNT OF OXYGEN IN CATHODE, G-MOLE/HR. 

ZH= REFORMER LENGTH, FT 
DINS(I): INLET AMOUNT OF GAS I, G-MOLE 
DNStl): INLET (OUTLET) AMOUNT OF GAS I, G-MOLE 
HA( J ) : SURFACE AREA OF HEAT EXCHANGER J 

HCAS(I),HCBS(I),HCCS(I): HEAT CAPACITY CONST. OF GAS I, CAL/G-MOLE 
OF THE FCRM:HCAS+HCBS*T+HCCS*T**2 
HS( I ) : HEAT OF FORMATION OF GAS I AT 293 K, 1 ATM 
NNS(I): STOICHIOMETRIC COEFFICIENT OF GAS I 
UMt I ) : MOLECULAR WEIGHT OF GAS I, G/G-MOLE 
DNSS! I , J ) : FLOW RATE OF GASJ IN STREAM I, G-MOLES/HR 

REAL KO , IDSH , IDTH 

DIMENSION GS( 7) , HS(7), HCAS(7>, HCBS(7), HCCS(7), WM(7) 

DIMENSION 0NSU7), DHSV(7), HLHVC7) 

DIMENSION DN5S( 39,7), T(39), P(39), TDNSSI 39 ) 

DIMENSION DNS1C7), DNS2(7), DHS(7), DHSM(7), DNSANt 7 ) , DNSCAI7), D 
1NSAI ( 7 ) , ONSC( 7 ) i DN5H ( 7 ) 

DIMENSION Q3T(10), HA(10), EFF(IO), NT(10) 

DIMENSION DNSR( 7 ) , DNSF(7) 

DIMENSION DPD( 2 ) > AHRN( £ ) , APPDI2), CLEPD(2), NTPDI2), SV(3) . 

DIMENSION AA1(2,7), AA2(2,7), AA3(4,7) 

NAMELIST /REPEN/ ZH , DX1 ,DX2 , DX3 ,KO , EA ,RHOB , EPS , S , DP ,DZZ 

NAMELIST /HE ATX/ CN,U,HA 

NAMELIST /OPFC/ TOPFC.UT.CD 

fJAMELIST /IHIT/ ONSM , TAT , PAT ■ SMRA, POPR 

NAMELIST /CON'D T/ IFUE L , ERR , IP , I , EXT , WAT , EXA 

NAMELIST /HEFDC/ NFH , NRH , BSPAC , ODTH , PITCH , CLH , IDSH , IDTH , FLOAR , SURF 
1C , CLENH , S1TS2 , DTH 

NAMELIST /PDSHH/ DPO , AHRN , APPO , CLEPD ,NTPD 

NAMELIST /PDFUH/ NTAF , FULE , WIDAF ,NPFU , NTAA , AIRL , WIDAA 

NAMELIST /TR/ TCPR 

NAMELIST /CATAI/ SRO,SA,CU,CL,ALFA,SN,FCONST,DKC 

COMMON /TC/ TC/CONS/A , B/GAG/GAGM , GAGA 
COMMON /CONST/ I 
COMMON /EXA/ EXA 

COMMON /ETHDA/ GS , HS ,HCAS ,HCBS , HCCS 
COMMON /HUMI/ WAT 
COMMON /EXT/ EXT 
COMMON /Ul/ U 

COMMON /CONFC/ E , ETH , El , EV, EC.EFC 
COMMON /HE/ HE 



0004360 
0004380 
0004400 
0004420 
0004440 
0004460 
0004430 
0004500 
0004520 
0004540 
0004560 
0004560 
0004600 C 
0004620 
0004640 C 
0004660 
0004680 
0004700 C 
0004720 
0004740 
0004760 C 
0004760 
00046C0 
0004820 C 
0004840 
0004360 
0004880 C 
C004900 
0004920 
0004940 C 
0004960 
C 004960 
0005000 C 
0005020 
0005040 
0005060 C 
0005060 
0005100 
0005120 C 
0005140 
0005160 
0005130 C 
0005200 
0005220 
0005240 C 
0005260 C 
0005280 C 
0005300 
‘0005320 
0005340 
0005360 
0005380 
0005400 
0005420 


COMMON /CN1/ CN 

COMMON /NAFH/ HNA,BPNA,WMNA,VHNA 
COMMON /WM/ UM 

COMMON /HEFOT/ NPH ,NRH ,BSPAC ,OOTH , PITCH .CLH , IDSH , IDTH , FLOAR ,SURFC , - 
1CLENH i S1TS2 i DTH 

COMMON /POSHT/ DFD , AHRN , APPO , CLEPO ,NTFO 

COMMON /PDFUT/ NTAF , NTAA , FULE , AIRL , WID AF .WIDAA.NPFU 

COMMON /REP/ KO ■ E A i RHOB , EPS i DZZ 

COMMON /SV/ SV.SVU 

COMMON /HLHV1/ HLHV 

COMMON /T'riCC/ AA1/VXPC/AA2/HTCPC/AA3 

COM.MON /CATAL/ SRO ,SA ,CU , CL , AIF A , SN , FCONST , ARE AF ,DKC 

DATA 0NSS/273*0./ 

READ THE OPERATION CONDITION OF FUEL CELL 
READ ( 5>0PFC ) 

WRITE ( 6 > OFFC ) 

READ THE INLET AMOUNT AND CONDITION 
READ (5.INIT) 

WRITE (6.INIT) 

READ OPERATION CONDITION 
READ (5.CCNDT) 

WRITE (6.C0NDT) 

READ OPERATING COEFFICIENT IN THE REFORMER! METHANE FUEL ONLY) 

IF (IFUEL.EO.l) READ (5, REPEN) 

IF (IFUEL.EQ.l) WRITE (6, REPEN) 

READ CONDITION OF HEAT EXCHANGER 
READ ( 5 >HEATX ) 

WRITE (6, HE ATX) 

READ. THE CONSTRUCTION OF HEAT EXCHANGER FOR CAL. PRESSURE DROP 
READ ( 5 >HE FDC ) 

WRITE ( 6 jHEFDC ) 

READ CONFIGURATION COEFF. OF SHIFT CONVERTER FOR CAL. PRESSURE DRO 
READ ( 5 > PDSHH ) 

WRITE (6, PDSHH) 

READ CONFIGURATION COEFF. OF FUEL CELL FOR CAL. PRESSURE DROP, 

READ ( 5 > PCFUH ) 

WRITE ( 6 , FDFUH ) 

READ THE OPERATING TEMP. OF REFORMER FOR FUEL NAPHTHA AND METHANOL 
IF ( IFUEL.NE . 1 ) READ (5,TR) 

IF (IFUEL.NE. 1) WRITE (6,TR) 

READ CATALYST CONSTANTS 
RE AO ( 5 i CATAI ) 

WRITE ( 6 iCATAI ) 


CHANGE THE THERMAL DATA FOR DIFFERENT FUEL INPUT 
IF (IFUEL.EQ.2) GO TO 1 
IF (IFUEL.EQ.3) GO TO 2 
GO TO 3 

1 GS( 1 ) = -38810 . 

HS(l)=-48050. 

HCAS( 1 )=4 , 394 
HCB5( 1 )-0 . 024274 
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0005440 

0005460 

0005400 

0005500 

0005520 

0005540 

0005560 

0005500 

0005600 

0005620 

0005640 

0005660 

0005600 

0005700 

0005720 

0005740 

0005763 

0005700 

0005800 

0005020 

0005040 

0005060 

0005800 

0005900 

0005920 

0005940 

0005960 

0005980 

0006000 

0006020 

0006040 

0006060 

0006080 

0006100 

0006120 

0006140 

0006160 

0006100 

0006200 

0006220 

0006240 

0006260 

0006280 

0006300 

0006320 

0006340 

0006360 

0006300 

0006400 

0006420 

0006440 

0006460 

0006400 

0006500 


HCCSf 1 )=-0 . 0000069 
WMt 11 = 32. 

HLHV( 1 )=-159258. 

AA3( 1 i 1 ) = 4 . 394 
AA3(2, 11=0. 013486 
AA3(3,l)=-2.1157E-06 
AA3( 4 i 1 ) =0 . 

AA1( 1 > 1 ) = 3. E-05 
AAl(2,l)=1.5E-02 
AA2( 1 >1 ) = 1 .8E-05 
AA2( 2 i 1 ) =0 . 0092 
GO TO 3 

2 GS( 11=6520. 

HS(l)=-42275. 

HCASt 11=7.094 
HCBSU 1=0.123447 
HCCS(l)=-0. 0000387 
WM ( 1 1 = 1 0 0 . 

HLHVt 1 1 = -1099580 . 

AA3(1, 11 = 7. 408 
AA3(2, 11=0. 062467 
AA3(3,l)=-1.0945E-05 
AA3(4>1 1=0. 

AA2(1, 11=2. 66 E-05 
AA2(2, 11=0. 0157 
AA1 (1,11=2. 5E-05 
AA1( 2 >1 1 = 0 . 7E-02 

3 CONTINUE 
IDNO=l 

C FUEL INPUT COMPRESSOR! PUMP)*** 

P( 1 ) = POPR*l . 02 

IF (IFUEL.EO.il CALL COMP ( DNSM, TAT, TOUT, PAT.Ptl) , POUM.GAGM, I, IP) 
IF (IFUEL.NE.l) CALL PUP ( DNSM ,TAT ,TOUT, PAT , P( 1 1 , POWN , I , IFUEL 1 
DO 301 IA=1 i I 

301 DNSSt 1 , IA 1 =DNSM( I A 1 

T(l)=TOUT 

C ASSUME THE COMP. OF 13TH FLOW 
DNSSt 13,21=0. 

DNSSt 13i 71=0 . 

IF (IFUEL. EQ. 2) GO TO 4 

IF (IFUEL. EO. 3) GO TO 5 

DNSSt 13)11=0. 147*DNSS( 1 j 1 1 

DNSS( 13,3 1 = 0 . 0018*DNSS( 1 , 1 )+DNSS( 1 ,3) 

DNSSt 13,41=0 . 85*DNSS( 1 , 1 )+ONSS( 1 ,4) 

DNSSt 13,5 1=3. 4*DNSS(1,1)+0NSS( 1,5) 

DNSSt 13,6 )=SMRA*DNSS( 1,11-1. 7*ONSS( 1,1) 

GO TO 6 

4 DNSSt 13,1 1=0. 00086 3487*DNSS( 1,1) 

DNSSt 13,31=0.01 

DNSSt 13,41=0 . 999128*DNSS( 1,1) 

DNSSt 13,5 1=2. 9974*DNSS( 1,1) 

DNSSt 13,61=0. 60087*DNSS( 1,1) 

GO TO 6 

5 DNSSt 13,1 1=0 . 



0006520 
0006540 
0006560 
0006580 
0006600 C 
0006620 
0006640 C 
0006660 
0006680 
0006700 
0006720 
0006740 
0006760 
0006760 
0006800 
0006820 C 
0006840 
0006360 
0006680 
0006900 
0006920 
0006940 
0006960 
0006980 
0007000 
0007020 
0007040 
0007060 
0007080 C 
0007100 C 
0007120 
0007140 
0007160 
0007180 
0007200 
0007220 
0007240 
0007260 
0007280 
0007300 
0007320 C 
0007340 
0007360 
0007360 
0007400 
0007420 
0007440 
0007460 
“0007480 
0007500 
0007520 C 
0007540 
0007560 
0007580 


DNSS( 13,31=0 . 0221*DNSS( 1 > 1 ) 

DNSSI 13)41=6. 971*DNSS( 1 > 1 ) 

DNSSl 13,51=21. 956*DNSS< 1,1) 

DNSSI 13,6 )=SMRA*DNSS( 1 , 1 1-13 . 9645*0NSS( 1,1 ) 

ASSUME THE PRESSURE OF 13TH FLOW 

6 P( 13 )=POPR*0 .958 

ASSUME THE PRESSURE OF 31ST FLOW 

7 IF (IDNO.GT.l) GO TO 8 
P( 31 ) = P( 13 1 

T( 31 )=TAT 

8 PINF=P( 13 1 
PINA=P( 31 1 
POPFC=P( 31 1 
DO 9 IB=1 , I 

9 0HS( IB 1 =DNSS( 13 1 IB 1 

FUEL CELL MASS BALANCE*** 

CALL FUCEI DNS , TOPFC , POPFC ,0NSAN ,DNSCA >DSO ,D5N ,DSHO >UT , I , PINF , PINA - 
l.IFUEL) 

P( 14 1 =PINF 

P( 32 1 =PINA 

DO 10 IC=1.I 

DNSSI 32 , IC )=DNSCA( IC 1 

10 DNSS( 14 , IC 1 =DHSAN( IC 1 
T( 32 1 =TOPFC 

T( 14 )=TOPFC 

IF ( IFUEL . EQ . 2 1 CK=1.5 

IF (IFUEL. EQ.l) CK=2 . 

IF (IFUEL. EQ. 3) CK=15. 

CAL. THE DIVIDER FACTOR FOR THE ASSUMPTION OF DEFINED EXTRA AIR IN 
NER 

GARM=1 ./( < ( 1 . +EXT*0 . 01 )*( ( DNSSI 14,3 )+DHSS( 14,511/2. +CK*DNSS( 14,1)1- 
1 l/DSO+l . 1 
DO 11 10=1,1 

11 DNSSI 28, ID 1=0. 

DNSSI 28,2) =DSO/GARM 
DNSSI 28,6 )=DSHO/GARt1 
DNSSI 28,7 )=DSN/GARM 
PI 28 )=PAT 

DO 12 IE=1 , I 

12 DNSAI( IE )=DHSS( 28 , IE 1 
AIR COMPRESSOR*** 

T( 28 ) = TAT 
POUT=POPFC*l .005 

CALL COMP I DNSAI ,T( 28 1 ,TOUT ,P( 28 1 , POUT , POWA ,GAGA,I , IP 1 
PI 29 ) = POUT 
DO 13 18=1,1 

13 DNSSI 29,18) =DNSAI ( 18 1 
DO 14 IG=1,I 

14 DNSI IG)=DNSS( 29, IG) 

TC 29 )=TOUT 
DIVIDER*** 

CALL DIVID (T(29),T(30),T(33) , DNS, DNSI, DNS2 ,GARM,I 1 

DO 15 IH=1,I 

DNSSI 33, IH)=DHS2(IH) 
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0007600 ONSSt 30, IH)=DNS1(IH) 

0007620 15 0NSS(31,IH)=0NSS(30,IH) 

0007640 P( 30 ) = P ( 29 ) 

0007660 P( 33 ) = P( 29 ) 

0007680 C CAL. THE 33RD FLOW COMP. UNDER THESE ASSUMPTION: 

0007700 C (1). COMPLETE COMBUSTION 

0007720 DNSS(33,2)=(1.+£XT*0.01)*< < ONSSt 14, 3 )+DNSS( 14,5 ) )/2 . +CK*DNSS( 14, 1 )- 

0007740 1) 

0007760 DNSS( 33,7 )=ONSS( 33,2)*3.76 

0007780 DNSS(33,6)=( < ONSSt 33 , 2 )+DNSS( 33, 7) )*28 . 6 )*WAT/18. 

0007800 T(33)=TAT 

0007820 DO 16 11=1,1 

0007840 DHS1 ( II )=ONSS( 14,11) 

0007360 16 DNS2 ( II ) =DNSS( 33,11) 

0007880 C MIXER*** 

0007900 CALL DMIX ( DNS1 ,DNS2 ,DNS , T( 14 ) , T< 33 ) , T( 15 ) , I , P( 14 ) , P( 33 ) , POUT ) 

0007920 P( 15 )=POUT 

0007940 DO 17 IJ=1,I 

0007960 DNSS( 15 , 1 J )=DNS( IJ ) 

0007980 17 DNSSt 16 , 1 J )=DNSS( 15 , IJ ) 

0008000 C CAL. THE COMP. OF 17TH FLOW UNDER THE SAME ASSUMPTION AS BEFORE 
0008020 IF (IFUEL.NE.3) DNSSt 17 ,4 )=0NSS( 15 , 3 )+DNSS( 15 ,4 ) +DNSS( 15 , 1 > 

0008040 IF (IFUEL.EQ.3) DNSSt 17 , 4 )=ONSS( 15 , 3 )+DNSS( 15 ,4 ) + 7 . *DNSS( 15, 1 ) 

0008060 IF (IFUEL.NE.3) ONSSt 17,6 )=DNSS( 15,5 )+ DNSSt 15,6 ) + 2 .*DNSS( 15 , 1 ) 

0003060 IF (IFUEL.EQ.3) DNSSt 17,6 ) =DNSS( 15 ,5 ) +DNSS( 15 , 6 ) +8 . *DNSS( 15 , 1 ) 

0008100 DNSSt 17,2 ) =ONSS( 15,2 )-0.5*( DNSSt 15, 3) + DNSSI 15,5 ))-CK*( DNSSt 15,1)) 

0008120 DNSSt 17,7) =DNSS( 15,7) 

0008140 00 18 IK=1,I 

0003160 ONSSt 18 , IK )=DNSS( 17, IK ) 

0008180 DNSSt 19 , IK ) =DNSS( 16 , IK ) 

0008200 DNSSt 20 , IK ) =DNSS( 19 , IK ) 

0008220 18 CONTINUE 

0008240 C CAL. THE COMP. OF 21ST FLOW 

0008260 DO 19 IL=1,I 

0008280 19 DNSSt 21 , IL )=DNSS( 20 , IL )+DNSS( 32 , IL ) 

0008300 TDNSSt 21 ) = 0. 

0008320 DO 20 IM=1,I 

0008340 20 TDNSSt 21 ) = TDNSSt 21 )+DNSS( 21 , IM ) 

0008360 C CAL. TEMP. OF 26TH FLOW 

0008380 C ASSUMPTION: 

0008400 C (1). SATURATED PRESSURE IS ESTIMATED BY (EXP(A-B/T)) 

0008420 C (2). UNDER SATURATED CONDITION 

0008440 C ASSUME THE PRESSURE OF 26TH FLOW IS IDENTICAL TO 25TH FLOW (THE PRESS. 
0008460 C DROP IS SMALL) 

0008480 P(25)=POPR*1.0005 

0008500 Pt 26 ) = P( 25 ) 

0008520 T( 26 )=B/( A-ALOG( Pt 26 ) ) ) 

0008540 DNSSt 26 ,6 ) =DNSS( 1 , 1 !*SMRA 

0008560 DO 21 IN=1,I 

0008530 21 DNSSt 25 , IN ) =DNSS( 26 , IN ) 

0008600 DO 22 10=1,1 

0008620 IF (IFUEL.NE.2) DNSSt 27,10 )=DNSS( 26 , 10 ) 

0008640 22 DNSSt 2 , 10 ) =ONSS( 1 , 10 ) 

0008660 DO 23 19=1,1 
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0008680 IF (IFUEL.NE.2) DNSS( 3 , 19 )=DNSS( 2 , 19 )+DNSS( 27,19 ) 

0008700 IF (IFUEL.EQ.2) DN5S( 3 , 19 ) =0NSS( 2 , 19 ) +DNSS( 26 , 19 ) 

0008720 IF (IFUEL.NE.2) DNSSt 4,19 )=DNSS( 3,19) 

0008740 IF (IFUEL.EQ.2) DNSS( 5 , 19 )=DHSS( 3 , 19 ) 

0006760 IF (IFUEL.NE.2) DHSS( 5 , 19 ) =DNSS( 4 , 19 ) 

0003780 23 CONTINUE 

0006800 C ASSUME THE TEMP. AND PRESSURE OF 5TH FLOW 

0008820 IF (IDNO.GT.I) GO TO 24 

0008840 IF (IFUEL.EQ.l) T(5)=1110.88 

0008860 IF (IFUEL.EQ.2) T(5)=661.59 

0008880 IF (IFUEL.EQ.3) T(5)=I052. 

0008900 P(5) = POFR 

0008920 24 CONTINUE 

0008940 IF (IDNO.GT.I) GO TO 25 

0008960 C ASSUME THE TEMP. OF 18TH FLOW 

0008960 IF (IFUEL.EQ.l) T(18)=1222. 

0C09000 IF (IFUEL.EQ.2) T(18)=863.88 

0009020 IF (IFUEL.EQ.3) T( 18 ) =1398 . 94 

0009040 C ASSUME P( 16 ) 

0009060 P< 16 ) =P( 15 )*0 . 99 

0009030 C DESIGN THE TUBE TO LET PRESSURE DROP THROUGH BURNER BE 6X 
0009100 6C4 F(17) = P(16)*(l.-0.06) 

0009120 C ASSUME PRESSURE DROP OF COMBUSTIAN GAS THROGH REFORMER TO BE 5 /. 

0009140 P( 18 ) = P( 17 )*0 . 95 

0009160 25 CONTINUE 

0009180 C E-3*** 

0009200 DO 26 IQ=1,I 

0009220 DNSHt IQ )=DNSS( 17, IQ ) 

0009240 26 DN3C( IQ )=DNSS( 5 , IQ ) 

0009260 PT=P( 5 ) 

0009280 PS=P( 18 ) 

0009300 IF (IFUEL.NE.2) CALL HEXC ( T( 18 ) ,DNSH ,DNSC,T( 4 ) ,T( 19 ) ,T( 5 ) ,QT,1 , HA- 

0009320 1( 3 ) , 3, 2 > I >PT, PS ,NT( 3 ) > IFUEL ) 

0009340 IF (IFUEL.EQ.2) CALL HEXC ( T( 18 ) ,DNSH ,DNSC ,T( 3 ) ,T( 19 ) , T( 5 ) ,QT , 1 , HA- 

0009360 1(3), 3, 2, I, PT, PS, NT (3), IFUEL) 

0009380 IF (IFUEL.NE.2) P(4)=PT 

0009400 IF (IFUEL. Eq. 2) P(3)=PT 

0009420 P( 19 ) = FS 

0009440 EFF( 3 )=HE 

0009460 QQT( 3 )=QT 

0009480 C E-4*** 

0009500 DO 27 IR=1,I 

0009520 DNSH ( IR ) =DNSS( 19 , IR ) 

0009540 DN5C( IR )=DNSS( 15,IR ) 

0009560 27 CCNTINUE 

0009580 PT=P( 15 ) 

0009600 PS=P( 19 ) 

0009620 CALL HEXC ( T( 19 ) ,DNSH ,DNSC ,T( 15 ) ,T( 20 ) , T( 16 ) ,QT , 1 ,HA( 4 ) ,4 , 1 , 1 , PT , P- 

0009640 IS, NT(4), IFUEL) 

0009660 I F ( ABS( ( PT-P( 16 ) )/( PT+P( 16 ) ) ) . LT. ERR ) GO TO 605 

0009680 P(16) = (PT+P(16) )/2. 

0009700 GO TO 604 

0009720 605 P(16)=PT 

0009740 P( 20 ) = PS 
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0009760 
0009780 
0009800 C 
0009820 
0009890 
0009860 
0009880 
0009900 
0009920 
0009990 C 
0009960 
0009980 
0010000 
0010020 
0010090 
0010060 
0010080 C 
0010100 
0010120 
0010190 
0010160 
0010160 
0010200 
0010220 
0010290 
0010260 C 
0010280 
0010300 
0010320 
0010390 
0010360 
0010380 
0010900 
0010920 C 
0010990 
0010960 
0010980 
0010500 
0010520 
0010590 
0010560 
0010580 
0010600 
0010620 
0010690 
0010660 
0010680 
0010700 C 
0010720 
0010790 
0010760 
0010780 
0010600 
0010820 C 


EFF( 9 )=HE 
QQT! 9 )=QT 

BUPNER ENERGY BALANCE*** 

DO 28 IS=1.I 

28 DNS(IS)=0NSS(16,IS> 

CALL FLAME ( DNS ,T( 16 ) , T( 17 ) , I , IFUEL ) 

DO 29 IT=1,I 

29 DNSS! 17 > IT ) =DNS( IT ) 

IF (IFUEL. NE.l) GO TO 33 
REFORMER! METHANE )*** 

NOR=DNSS( 1,1 1/60. 

DO 30 IU=1,I 

DNSR! IU)=DNSS( 5, IU 1/953. 6/NOR 

30 ONSF ( IU 1 =DNSS( 17 , IU 1/953 . 6/NOR 

CALL KREF ( DNSR >0NSF>0X1,DX2,DX3,P(5)>T(5)>T(17)>ZH,P(7),T(7) 

1 , TTEST ,S,DP,1 1 

TEST THE ASSUMPTION OF 18TH FLOW 

IF ( ABS! ( TTEST-T! 18 ) )/( TTEST+T! 18 ) ) ) . LT . ERR ) GO TO 31 

T(181=TTEST 

GO TO 25 

31 DO 32 IV=1,I 

DNSS! 7 , IV 1 =DNSR( IV 1*953 . 6*N0R 

32 DNSS! 18 , IV )=DNSF ( IV 1*953 . 6*N0R 
T(18)=TTEST 

GO TO 39 

E-8! METHANOL AND NAPHTHA)*** 

33 DO 39 IW=1,I 

DNSH! IW )=DNSS( 17 , IW ) 

39 DNSC! IW 1 =DNSS( 5 , IW 1 
PTEST=P( 17 1 
PS = P( 5 1 

CALL HEXC ( T( 17 1 ,ONSH ,DNSC ,T( 5 1 , TTEST ,T( 6 ) ,(?T , 2 ,HA( 8 1 , 8, 1 , 1 , PTEST , 
1FS,NT(8), IFUEL) 

TEST THE ASSUMPTION OF 18TH FLOW 

IF ( ( ABS! ( TTEST-T! 18 1 )/( TTEST+T! 18 )) ) . LT. ERR 1 .AND . ( ABS! ( PTEST-P! 18 
1 1 )/( PTEST+P! 181)1. LT.ERR 1 1 GO TO 35 
T(18)=(TTEST+T(18))/2. 

P( 18 ) = ( P( 18 l+PTEST 1/2 . 

GO TO 25 

35 QQT(8)=0T 
T( 18 )=TTEST 
PI 18 )=PTEST 
PS 6 ) = PS 
EFF(8)=HE 

DO 36 IX=1,I 
DNSS! 6 , IX l-ONSC! IX 1 

36 DNSS! 18,IX)=DNSH< IX) 

REFORMER! METHANOL AND NAPHTHA)*** 

DO 37 I Y=1 , I 

37 DNS! IY ) -DNSS! 6 , 1 Y 1 

CALL ENRE ( DNS , TOPR , POPR ,T( 6 ) , T( 7 ) , I , IP, IFUE L ) 

DO 38 12=1,1 

38 DNSS! 7,IZ)=0NS(IZ) 

CAL. PRESSURE DROP 



0010640 
0010660 
0010660 
0010900 
0010920 
0010940 C 
0010960 
0010980 
0011000 
0011020 
0011040 
0011060 
0011080 
0011100 
0011120 
0011140 
0011160 
0011160 
0011200 
0011220 
0011240 
0011260 
0011280 C 
0011300 
0011320 
0011340 
0011360 
0011380 
0011400 
0011420 
0011440 
0011460 
0011480 
0011500 
0011520 
0011540 
0011560 
0011580 
0011600 
0011620 
0011640 
0011660 
0011680 C 
0011700 
0011720 
0011740 
0011760 
0011780 
'0011800 
0011820 
0011840 
0011860 
0011880 
0011900 


CALL PDSH (0NS,P(6),P(7),T0PR,2,IFUEL) 

IF ( IFUEL. Eq. 2 ) GO TO 56 

39 DO 40 IAA=1,I 

DNSH( IAA)=DNSS(7,IAA) 

40 DNSC(IAA)=DNSS(3,IAA) 

E-2*** 

PT=P( 4 ) 

PS-P( 7 ) 

IF (IFUEL. Eq.l) CALL HEXC(T(7),0NSH,DNSC,T(3),T(8),T(4),qT,l,HA(2),- 
12 > 2 > I . PT i FS ,NT( 2 ) > IFUEL ) 

IF( IFUEL . Eq. 3 ) CALL HEXC( T( 7 ) ,0NSH ,DNSC ,T( 3 ) .T( 9 ) ,T( 4 ) , qT , 1 .HANK- 
IE > 2 > I > PT i FS > NT( 2 ) > I FUEL ) 

P( 3 )=PT 

IF (IFUEL. EQ.l) P( 8 ) = PS 
IF (IFUEL. EQ. 3) P(9)=PS 
EFF( 2 )=HE 
qqT(2)=qT 
DO 41 IBB=1,I 

IF (IFUEL. Eq.l) DNSS(8,IBB)=DNSH( IBB) 

IF (IFUEL. Eq. 3) DNSS(9,IBB)=0NSH(IBB) 

41 DNSS! 4 i IBB ) =DNSC( IBB ) 

IF (IFUEL. Eq. 3) GO TO 46 
HIGH TEMP. SHIFT CONVERTER*** 

T0FHS=T(8) 

42 CONTINUE 
F0FHS=P<8) 

DO 43 ICC=1,I 

43 DNS( ICC )=DNSS( 8i ICC ) 

CALL ENSH ( DNS ,T( 8 ) ,T( 9 ) ,TOPHS , POPHS , I , IP , IFUEL ) 

P(9) = F0PHS 
TM=(T(9)+T(8) )/2. 

IF ( ( IP. Eq . 2 ) . OR . ( ABS( ( TM-TOPHS )/( TM+TOPHS ) ) . LT . ERR ) ) GO TO 44 

TOPHS=TM 

GO TO 42 

44 CONTINUE 

DO 45 IDD=1,I 

45 DHSS(9,IOD)=ONS(IDD) 

46 CONTINUE 

DO 47 IEE=1,I 
DNSH(IEE)=DNSS(9,IEE) 

47 DNSC( I EE )=DNSS( 1 >IEE ) 

IF (IFUEL. Eq. 3) GO TO 48 
E-l( METHANE AND NAPHTHA)*** 

PT=P(1) 

PS=P( 9) 

CALL HEXC (T(9),DNSH,DNSC,T(l),THO,TCO,qT,l.HA(l),l.l,I,PT,PS,NT(l- 
1), IFUEL) 

P(10)=PS 
P( 2 ) = PT 
EFF( 1 )=HE 
qqT( 1 )=QT 
T( 2 )=TCO 
T( 10 )=THQ 
GO TO 52 
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0011920 98 CONTINUE 

0011990 C DESIGN THE HEAT EXCHANGER 1 TO VAPORIZE NAPHTHA AND RISE THE TEMP. 

0011960 C TO 900 

0011980 C K FOR NAPHTHA FUEL 

0012000 T( 2)=900. 

0012020 QQT( 1 ) =HNA*252 . /9S3 . 6*100 . *1 .8*DNSS( 1 > 1 )*( BPNA+273 .—TCI)) +DNSS( 1,1- 

0012090 1 )*VHNA+HCASt 1 )*< T< 2 )-T( 1 ) )+HCBS( 1 )*( T< 2 )**2-T( 1 )**2 >/2 .+HCCS( 1 )*( T- 

0012060 2( 2 )**3-T( 1 )##3 )/3. 

0012080 C ASSUME EFFICIENCY OF H-E 1 IS 0.7 
0012100 EFF( 11=0.7 

0012120 HA(1)=CN*QQT(1)/EFFC1)/(T(9)-T(1))/U 

0012190 C ASSUME AVERAGE TEMP. OF HOT SIDE 

0012160 TAVG=T(9)-50. 

0012180 99 CHH=0. 

0012200 00 50 IB=1 > I 

0012220 50 CHH=CHH+DNSS( 9, IB )#( HCAS( IB )+HCBS( IB )*TAVG+HCCS( IB }*TAVG*#2 ) 

0012290 T( 10 )=T( 9 )-GQT( 1 l/CHH 

0012260 C TEST THE ASSUMPTION OF AVERAGE TEMP. 

0012280 IF ( ABS( <T(10)+T(9))/2. -TAVG ) . LT . ABS( (T(10)+T(9))/2. +TAVG 1*0 ■ 001 J - 

0012300 1GO TO 51 

0C12320 TAVG= ( T( 9 ) +T( 10 ) 1/2 . 

0012390 GO TO 99 

0012360 51 CONTINUE 

0012380 C ASSUME PRESSURE DROP IS NEGLECTABLE 

0012900 P< 2 )=P( 1 ) 

0012920 P( 10 ) = P( 9 ) 

0012990 52 CO 53 IFF=1,I 

0012960 IF ( IFUEL. EO. 3 ) DNSS( 10 , IFF ) =DNSS( 9 , IFF ) 

0012980 53 IF (IFUEL.EQ.l) DNSS( 10 , IFF )=DNSH( IFF ) 

0012500 DO 59 IGG=1,I 

0012520 DNSH< IGG )=DHSS( 10 , IGG ) 

0012590 59 DNSC< IGG )=DNSS( 26 > IGG ) 

0012560 C E-9*** 

0012580 PT=P( 26 ) 

0012600 PS=P(10) 

0012620 CALL HEXC ( T( 10 ) ,DNSH,DNSC,T( 26 ) ,T( 11 ) ,T( 27 ) ,qT, 1 ,HA( 9 ) , 9, 1 , I ,PT, P- 

0012690 IS , NT( 9 ) , IFUEL ) 

0012660 P(11) = PS 

0012680 P( 27 )=PT 

0012700 EFF(9)=HE 

0012720 QGT( 9 )=QT 

0012790 DO 55 IHH=1,I 

0012760 DNSS(11,IHH)=DNSH(IHH) 

0012780 DNSS(27,IHH)=DNSC(IHH) 

0012800 DNS1 ( IHH ) =DNSS( 2 i IHH ) 

0012820 DNS2 ( IHH ) = DNSS( 27 , IHH ) 

0012890 55 CONTINUE 

0012860 GO TO 60 

0012880 C E-l( METHANOL )**» 

0012900 56 DO 57 111=1,1 

0012920 DNSH (III )=DNS5( 7, III ) 

0012990 57 DNSC(III)=DNSS(1,III) 

0012960 PT=P( 1 ) 

0012980 PS=P<7) 



0013000 
0013020 
0013040 
0013060 
0013030 
0013100 
0013120 
0013140 
0013160 
0013130 
0013200 
0013220 
0013240 
0013260 
0013230 C 
0013300 
0013320 
0013340 
0013360 C 
0013330 
0013400 
0013420 
0013440 C 
0013460 
0013460 
0013500 
0013520 
0013540 
0013560 
0013560 
0013600 
0013620 
0013640 
0013660 
0013680 
0013700 
0013720 
0013740 
0013760 
0013780 
0013800 
0013820 C 
0013840 
0013860 
0013880 
0013900 
0013920 
0013940 
'0013960 
0013980 
0014000 
0014020 
0014040 
0014060 


CALL HEXC (T(7),DNSH,DNSC,T(l),TH0,TC0,QT,l,HA(l),l,l,I,PT,PS,m'(l- 
1), IFUEL) 

PTE=PS 
P(2) = PT 
EFF( 1 )=HE 
QQT( 1 )=QT 
DO 58 IJJ=1,I 

58 DNSS( 2 i XJJ )=DNSC( IJJ ) 

T(2)=TCO 

T( 13 )=THO 

DO 59 IJJ=1,I 

DUSK IJJ )=0NSS( 2 >IJJ ) 

DHS2( IJJ )=DNSS( 26 , IJJ 1 

59 CONTINUE 
MIXER*** 

CALL DMIX (DNS1,0NS2,DNS,T(2),T(26),TESTT,I,P(2),P(26),PTEST) 

GO TO 61 

60 CALL DMIX ( DNS1 , DNS2 , DNS ,T< 2 ! ,T( 27 ) , TESTT , I , P( 2 ) , P( 27 ) , PTEST ) 

TEST TEMP. OF 3RD FLOW 

61 IF ( ( ABS( ( TESTT-T( 3 ) )/( TESTT+Tt 3 ) ) ) . LT.ERR ) . AND . ( ABS( ( PTEST-Pl 3 ) )/- 
1 ( PTEST +Pt 3 ) ) ) . LT . ERR ) ) GO TO 70 

IF ( IFUEL . EQ . 2 ) GO TO 66 

DO 62 IKK=1 t I 
DMSS( 3, IKK )=DNS( IKK ) 

DNSH ( IKK )=DNSS( 7 . IKK ) 

62 DNSCC IKK ) r DNSS( 3»IKK ) 

P(3) = (P(3) + PTEST)/2. 

T( 3) = ( TESTT + T( 3 I )/2 . 

PT=P( 3) 

PS=P(7) 

IF (IFUEL. EQ. 3) GO TO 63 

CALL HEXC (T(7),DNSH,DNSC,T(3),T(8),T(4),QT,1>HA(2),2,1,I,PT,PS,NT- 
1(2), IFUEL) 

P(4) = PT 
P(8)=PS 
GO TO 64 

63 CALL HEXC ( T< 7 ) , ONSH , DNSC ,T( 3 ) ,T( 9 ) ,T( 4 ) ,QT , 1 , HA( 2 ) , 2 , 1 , I , PT , PS , NT- 
1(2 ), IFUEL) 

P(4) = PT 
P< 9 )=PS 
E-3*** 

64 DO 65 ILL=1,I 

0NSS( 4 , ILL )=DNSC( ILL) 

DNSH (ILL )=DNSS( 18 > ILL ) 

65 DNSC( ILL ) -DNSS( 4 ■ I LL ) 

PT=P(4) 

PS=P(18) 

CALL HEXC (T(18),0NSH,DNSC.T(4),T(19),T(5),QT,1,HA(3I,3,1,I,PT,PS,- 
1NT(3), IFUEL) 

P( 5 )=PT 
P( 19 )=PS 
GO TO 68 

66 DO 67 IMM=1,I 
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0014000 

0014100 

0014120 

0014140 

0014160 

0014180 

0014200 

0014220 

0014240 

0014260 

0014200 

0014300 

0014320 

0014340 

0014360 

0014380 

0014400 

0014420 

0014440 

0014460 

0014480 

0014500 C 

0014520 

0014540 

0014560 

0014580 

0014600 

0014620 

0014640 C 

0014660 

0014680 

0014700 

0014720 

0014740 

0014760 

0014760 

0014800 

0014820 C 

0014840 

0014860 

0014860 

0014900 

0014920 

0014940 

0014960 

0014980 

0015000 

0015020 

0015040 

0015060 C 

0015080 

0015100 

0015120 

0015140 


ONSSt 3,IMM)=0NS( IMM) 

DNSH ( IMM ) =DNSS( 18 , IMM ) 

67 0NSC( IMM )=DNSS( 3 , IMM ) 

T( 3 ) = ( T( 3 HTESTT )/2 . 

PT-Pl 18) 

PS=( P( 3 J + PTEST )/2 . 

CALL HEXC (T(18),DNSH,DNSC,T(3),T(19),T(5),qT,l,HA(3),3,l,I,PT,PS,- 
1NT( 3 ) , IFUEL ) 

P ( 1 9 ) = PT 
P(5) = PS 
EFF( 3)=HE 
qqT( 3)-QT 

68 DO 69 INN=1,I 

69 DNSS< 5 > INN ) =DNSC( INN ) 

IDNO-IDNO-fl 

GO TO 24 

70 CONTINUE 

IF (IFUEL. Eq. 2) GO TO 77 

DO 71 IPP=1,I 

DNSHt IPP)=DNSS( 11 ,IPP ) 

71 DNSC(IPP)=DNSS(30,IPP) 

E-6*** 

PT=P( 30 ) 

PS=P(11) 

CALL HEXC <T(ll),0NSH,DNSC,T(30),T(12),T(31),qT,l,HA(6),6,l,I,PT,P- 
1S,NT(6), IFUEL) 

P( 12 ) = PS 
PTEST=PT 

TEST THE PRESSURE OF 31ST FLOW 

IF ( ABS( ( PTEST-P( 31 ) )/( PTEST+PC 31 ) ) ) . LT. ERR ) GO TO 72 
P C 31 ) = ( PTEST+P( 31 ) )/2 . 

GO TO 8 

72 EFF( 6 ) =H E 
QQT( 6 )=QT 

DO 73 IQQ=1,I 

ONSS( 12 , iqq ) =DNSH( iqq ) 

73 DNSS ( 31 1 IQq )=DNSC( IQQ ) 

LOW TEMP. SHIFT CONVERTER*** 

TOPLS : T( 12 ) 

POPLS=P( 12 ) 

74 CONTINUE 

DO 75 IRR=1,I 

75 DNS(IRR)=DNSS(12,IRR) 

CALL ENSH ( DNS , T( 12 ) ,T( 13) ,TOPLS , POPLS, I , IP, IFUEL ) 
TM=(T(12)+T(13))/2. 

IF ( ( IP . EO. 2 ) . OR . ( ABS( ( TM-TOPLS )/( TM+TOPLS ) ) . LT . ERR ) ) GO TO 76 

TOPLS=TM 

GO TO 74 

76 CONTINUE 

TEST THE ASSUMPTION OF 13TH FLOW 

77 DO 78 ISS=1,I 

IF ( (IFUEL. NE. 2). AND. (DNS(ISS).LT. 0.50). AND. (DNSSt 13, ISS).LT. 0.50)- 
1) GO TO 78 

IF ( ( IFUEL. EG. 2 ) .AND . (0NSH( IS5 ) . LT.0.50 ) .AND . (DNSSf 13,ISS).LT.0.50- 



0015160 

00151S0 

0015200 

0015220 

0015240 

0015260 

0015280 

0015300 

0015320 

0015340 

0015360 

0015380 

0015400 

0015420 

0015440 

0015460 

0015480 

0015500 

0015520 

0015540 

0015560 

0015580 

0015600 

0015620 

0015640 

0015660 

0015680 

0015700 

0015720 

0015740 

0015760 

0015780 

0015800 

0015820 

0015840 

0015860 

0015380 

0015500 

0015920 

0015940 

0015960 

0015980 

0016000 

0016020 

0016040 

0016060 

0016080 

0016100 

"0016120 

0016140 

0016160 

0016130 

0016200 

0016220 


1 )) GO TO 78 

IF ( ( IFUE L . NE . 2 ) . AND . ( ABS( ( DNS( ISS )-DNSS( 13 i ISS ) )/( DNS( ISS ) + DNSSt 1 
13 1 ISS ) ) ) . LT . ERR ) 1 GO TO 78 

IF ((IFUEL.EQ.2) .AND . ( ABS( ( ONSH ( ISS ) -ONSS( 13 , ISS ) )/( DNSH ( ISS ) +DNSS 
1(13, ISS))). LT. ERR)) GO TO 78 
GO TO 79 

78 CONTINUE 

IF ( IFUEL . EG. 2 ) POPLS=PTE 

IF ( ABS( ( POPLS-P< 13 ) )/( POPLS+Pt 13 ) ) ) . LT . ERR ) GO TO 81 

79 DO 80 ITT=1 , I 

IF (IFUEL. Ed. 2) DNSSt 13 , ITT ) = ( DNSS! 13 , ITT)+DNSH( ITT ) )/2 . 

80 IF (IFUEL. NE. 2) DNSSt 13 , ITT ) = ( DNSSt 13 , ITT HDNSt ITT ) )/2 . 

P( 13 ) = ( POPLS+Pt 13 ) )/2 . 

GO TO 7 

81 CONTINUE 

DO 82 IUU=1,I 

IF (IFUEL. Ed. 2) DNSSt 13, IUU )=DNSH( IUU ) 

82 IF (IFUEL. NE. 2) DNSSt 13 , IUU )=DNS( IUU ) 

TDN13=0 . 

DO 83 IVV=1,I 

83 TDN13=TDN13+DNSS< 13.IW) 

FCO=ONSS( 13,3 )/TDN13 

DO 84 IWW=1,I 

DNSANt IWN ) =0NSS< 13 , INK ) 

84 DNSCA( IWW )=DNSS( 31 > IWW ) 

FUEL CELL ENERGY BALANCE*** 

AREAF=FULE*AIRL*30. 48**2 
CD-CO/'IOOO . 

OU=l . /( EXA*0 . 01+1 . ) 

CALL ENFU( DN5AN , DNSCA,T( 31 ) , T( 13 ) ,TOUT ,TOPFC,POPFC > VOP ,UT , 1 , 1 ,QY , 
INK ,OU,CD ,CL > IFUE L ) 

QRT=QY 

CALCULATE THE OUTPUT OF AC POWER 

AC- ( -1 .0148+SQRTt 1 . 0148**2-4 . *0 . 0456/108.*! 0 . 0472*108 . -WK ) ) ) 

l/( 2. *0.0456/108. ) 

DO 86 IYY=1 , I 

86 DNSC( I YY )=DNSS< 25 , IYY ) 

DO 87 IZZ=1,I 

DUSK IZZ ) =DNSS( 20 , IZZ ) 

87 DNS2 ( IZZ )=DNSS( 32 ,IZZ ) 

MIXER*** 

CALL DMIX ( DNS1 ,DNS2 ,ONS,T( 20 ) ,T( 32 ) ,TOUT ,I,P(20),P(32), POUT ) 

P(21)=POUT 

DO 83 11=1,1 

DNSSt 21,11) =0NS( II ) 

88 DNSH( II )=0NSS( 21,11) 

T( 21 )=TOUT 

ASSUME P( 22 ) 

P( 22 ) = P( 21 )*1 .001 

602 T( 22 ) = -B/( ALOGt ( DNSSt 1 , 1 )*SMRA-DNSS( 21,6) )/(DNSS( 1 , 1 )*SMRA- 
1TDNSS(21))*P(22))-A) 

E-10*** 

T( 25 )=T( 22 ) 

CALL CD PH (QRT,T( 25 ) ,T( 26 ) ,DNSC,P( 25 ) ,1 ,QGT( 10 ) ) 
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0016240 
0016260 C 
0016280 
0016300 
0016320 
0016340 
0016360 
0016380 
0016400 C 
0016420 
0016440 
0016460 C 
0016480 
0016500 
0016520 
0016540 
0016560 
0016580 
0016600 
0016620 C 
0016640 
0016660 
0016660 
0016700 
0016720 
0016740 
0016760 
0016780 
0016800 
0016820 
0016840 
0016860 
0016880 C 
0016900 
0016920 
0016940 
0016960 
0016980 
0017000 
0017020 
0017040 C 
0017060 
0017080 
0017100 
0017120 
0017140 
0017160 
0017180 
0017200 
0017220 
0017240 C 
0017260 
0017280 
0017300 


GQT( 7 )=QRT 

ASSUME HOT WATER IS FEEOING IN AT 3330K AND FEEDING OUT AT 355 K 
DHSSl 34,6 )=QQT( 7)/< 355.-333. )/l ./18. 

DNSS(35,6)=0N3S(34,6) 

T( 341 = 333. 

T( 35 1 = 355 . 

P( 34 ) = PAT 
P( 35 )=PAT 
E-5*** 

CALL COND (T(211,T(22) ,DNSH >QT ,HCH ,1 1 
G3T( 5 )=QT 

ASSUME HOT WATER IS FEEOING IN AT TAT AND FEEDING OUT AT 3550K 
DTL=( (T( 21 1-355. )-( T( 22 1 -TAT 1 l/ALOG( ( T( 21 1-355 . )/( T( 22 l-TAT 1 1 
T( 36 )=TAT 
T( 37 1 = 355 . 

P( 361 = PAT 
P( 37 1 =PAT 

ONSS< 36,6 )=QGT( 51/1. /18./( 355. -TATI 
DNSS( 37,6) =DNSS( 36,6) 

ASSUME E-5 IS THE TYPE OF COUNTERFLOW 
HA( 5 )=QQT( 5 l/U/l ./DTL 
THM=(T( 36 l+T (37)1/2. 

TCM= ( T( 22 1 +T( 21 ) 1/2 . 

CALL HEPD( DNS , DNSH , THM , TCM ,HA( 5),P(21),P(36),5,DPJ,DP,NT,T(22), 
1T( 211,4) 

P22TE=P(21)-DPJ 

I F ( A3S( (P22TE-P(22))/( P22TE + P( 22)1). LT. ERR 1 GO TO 603 
P(22) = (P(22)*P22TE)/2. 

GO TO 602 
603 DO 89 12=1,1 

DNSSt 22,12 )=DNSH( 12 1 

89 DNSC 12 1 =DMSS( 22,12) 

SEPARATER*** 

POPS=P( 22) 

CALL SEPAR ( T( 22 1 ,POPS,T( 23 1 ,T( 24 1 .DNS.DNSL ,DNSV, 1 1 

DO 90 13=1,1 

DNSSt 24, I3)=DNSL(I3) 

90 DNSS( 23,13) =DNSV( 13 1 
P( 24 ) = POPS 

P< 23 ) = POPS 
PUMP*** 

DO 901 14=1,1 
901 DNS( 14 )=DNSS( 24,14) 

CALL FUMP ( ONS,T( 24 1 ,T( 25 1 , P( 24 ) ,P( 25 1 , POWS,I 1 
DO 92 IB=1 , 39 
TDNSS( IB 1 = 0. 

DO 91 IA=1 , I 

TDNSSt IB )=TDNSS( IB )+0NSS( IB.IA) 

91 CONTINUE 

92 CONTINUE 

WRITE THE RESULT 

IF ( IFUEL. EQ.l 1 WRITE (6,108) SMRA 
IF ( IFUEL. EQ. 2 ) WRITE (6,109) SMRA 
IF (IFUEL. Eq. 3) WRITE (6,110) SMRA 



0017320 

0017340 

0017360 

0017380 

0017400 

0017420 

0017440 

0017460 

0017480 

0017500 

0017520 

0017540 

0017560 

0017580 

0017600 

0017620 

0017640 

0017660 

0017680 

0017700 

0017720 

0017740 

0017760 

0017780 

0017800 

0017820 

0017840 

0017860 

0017680 

0017900 

0017920 

0017940 

0017960 

0017980 

0018000 

0018020 

0018040 

0018060 

0018080 

0018100 

0018120 

0018140 

0018160 

0018180 

0018200 

0018220 

0018240 

0018260 

'0018280 

0018300 

0018320 

0018340 

0018360 

0018380 


WRITE (6,103) P( 5 ) , P( 7 ) ,T( 5 ) , T( 7 ) 

IF ( IFUEL. EQ.l) WRITE (6,104) TOPHS, P( 8 ) ,T( 8 ) ,T( 9 ) 

IF (IFUEL. NE. 2) WRITE (6,105) TOPLS, P( 13 ) ,T( 12 ) ,T( 13 ) 

WRITE (6,106) T( 22 ) ,POPS 
WRITE (6,99) FCO 

WRITE (6,118) TOPFC , POPFC , E , VOP, CD >CL»UT , OU 
WRITE (6,102) T( 13 ) »T( 31 ) 

WRITE (6,100) ETH >EI > EV , EC 
WRITE (6,107) EFC 
WRITE (6,101) WK,QRT 
WRITE (6,117) AC 

C WRITE THE MATERIAL IN JTH FLOW 
18280 WRITE (6,119) 

IF (IFUEL. EQ.l) WRITE (6,94) 

IF (IFUEL. EQ. 2) WRITE (6,95) 

IF (IFUEL. EQ. 3) WRITE (6,96) 

IF (IFUEL. EQ.l) WRITE (6,93) ( ( ( J , ( DNSS( J ,IA ) , IA=1 ,1 ) ,TDNSS( J ) ,T( J 

1),P( J)),J=1,5),((II,(DNSS(II,IA),IA=1,I),TDNSS(II),T(II),P(II)),II 
2 = 7,37) ) 

IF (IFUEL. EQ. 2) WRITE (6,93) ( ( ( J , ( DNSS( J , IA ) , IA=1 , 1 ) , TDNSSt J ) , T( J 
1),P< J)),J=1,3),((II,(DNSS<II,IA),IA=1,I),TDNSS(II),T(II),P(II) ),II 
2=5,7) ) 

IF (IFUEL. EQ. 3) WRITE (6,93) ((( J ,( DNSS( J , IA ), IA=1 , 1 ), TDNSSt J ) ,T( J 
1),P(J)),J=1,7),((II,( DNSS( II,IA),IA=1,I) , TDNSSt II),T(II),P(II)),II 
2=9,37) ) 

C WRITE DUTY OF HEAT EXCHANGER OR CONDENSER 
WRITEt 6 ,119 ) 

WRITE (6,111) 

WRITE (6,97) 

WRITE (6,98) (QQT( IA) ,IA=1 >10 ) 

C WRITE SURFACE AREA OF HEAT EXCHANGER 
WRITE (6,112) 

WRITE (6,97) 

WRITE (6,98) ( H A( IA ) , IA=1 , 10 ) 

C WRITE THE EFFICIENCY OF HEAT EXCHANGER 
WRITE (6,113) 

WRITE (6,97) 

WRITE (6,98) (EFF(IA),IA=1,10) 

IF (IFUEL. EQ.l) WRITE (6,114) POWA , POWM , ROWS 
IF (IFUEL. EQ. 2) WRITE (6,115) POWA, POWN, PONS 
IF (IFUEL. EQ. 3) WRITE (6,116) POWA , FOWN, PCWS 
STOP 
C 


93 FORMAT ( IX, 2X , 12 , 2X, F8. 2 ,1X, F8. 2 , IX, 3X, F8. 2 ,6X, F8. 2 ,5X , F8. 2 , F10 . 2 , 
1 IX , F8 . 2 , 3X , F 10 . 2 , 4X , F7 . 2 , 2X , F6 . 4 ) 

94 FORMAT (IX, 'STREAM METHANE OXYGEN CAR. 

1H YDROGEN WATER NITROGEN FLOW RATE 

95 FORMAT (IX,' STREAM METHANOL OXYGEN CAR. 

1HYDROGEN WATER NITROGEN FLOW RATE 

96 FORMAT ( IX,' STREAM NAPHTHA OXYGEN CAR. 


1HYDRCGEN WATER 
97 FORMAT (IX,' 

1 E-5 E- 

2 E-10 ’ ) 


NITROGEN 
E- 1 


FLOW RATE 
E-2 

E-7 


MONOXIDE 
TEMP.(K) 
MONOXIDE 
TEMP.(K) 
MONOXIDE 
TEMP.(K) 
E-3 

E-8 


CAR. DIOXIDE 
PRE.(ATM) ' ) 
CAR. DIOXIDE 
(ATM)' ) 
DIOXIDE 
(ATM) ' ) 
E-4 


FRE. 

CAR. 

PRE. 


E-9 
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0018400 

0018420 

0018440 

0016460 

0018480 

0018500 

0018520 

0018540 

0013560 

0018530 

0018600 

0018620 

0018640 

0018660 

0016680 

0018700 

0018720 

0018740 

0016760 

0018780 

0018600 

0018820 

0018340 

0018860 

0018880 

0018900 

0018920 

0018940 

0018960 

0016980 

0019000 

0019020 

0019040 

0019060 

0019080 

0019100 

0019120 

0019140 

0019160 

0019180 

0019200 

0019220 

0019240 

0019260 

0019280 

0019300 

0019320 

0019340 

0019360 

0019380 

0019400 

0019800 

0019900 

0020000 


98 FORMAT ( IX, IOC E12 .5, IX )/ ) 

99 FORMAT ( IX. 'THE FRACTION OF CO IN THE FEED IS',F8.5/) 

100 FORMAT (IX, 'THE THERMODYNAMIC EFFICIENCY OF FUEL CELL IS' ,E13.5, 'T- 
1HE CURRENT EFFICIENCY IS',E13.5/' THE VOLTAGE EFFICIENCY IS',E13.5- 
2 ■ 'THE HEATING VALUE EFFICIENCY IS',E13.5) 

101 FORMAT (IX, 'THE ELECTRICAL WORK IS' ,E13.5, 'KW V THE TOTAL HEAT R- 
1ELEASE IS' ,E13.5, 'CAL'//) 

102 FORMAT (IX, 'THE ANODE SIDE INLET TEMP. IS',F7.2,' K'/' THE CATHODE- 
1 SIDE INLET TEMP. IS',F7.2,' K') 

103 FORMAT C THE REFORMER IS OPERATING UNDER THESE CONDITIONS ',6X,/'- 

1 INLET PRESSURE' ,F7. 2, 'ATM' , ' OUTLET PRESSURE ' ,F7. 2, 'ATM'/ ' - 

2 INLET TEMP. :',F7.2,' K',' OUTLET TEMP. :',F7.2,' K'/l 

104 FORMAT ( ' THE HIGH TEMP. SHIFT CONVERTER IS OPERATING UNDER THESE - 
1CONDITIONS ' >6X,/ 

2 ' OPERATING TEMP. : ' , F7 . 2 , ' KV OPERATING PRESSURE :', F7. 2 , 'ATM- 

3'/' INLET TEMP. : ' ,F7.2, ' K'/' OUTLET TEMP. : ' ,F7. 2 , ' KV) 

105 FORMAT (IX, 'THE LOW TEMP. SHIFT CONVERTER IS OPERATING UNDER THESE- 

1 CONDITIONS ' ,6X,/ 

2 ' OPERATING TEMP. : ' ,F7. 2 , ' KV OPERATING FRESSURE : ' ,F7.2, 'ATM- 
3V INLET TEMP. : ' ,F7.2, ' K'/' OUTLET TEMP. : ' , F7. 2 , ■ KV) 

106 FORMAT ( ' THE LIQUID SEPERATER IS OPERATING UNDER THESE CONDITIONS- 
1 1 ,6X,/ 1 OPERATING TEMP . : ' , F7. 2 , * KV OPERATING FRESSURE :', F7. 2 ,' - 
2ATMV) 

107 FORMAT (IX, 'THE FUEL CELL EFFICIENCY IS',F6.4) 

108 FORMAT Cl', IX, 'THE STEAM/METHANE RATIO IN THE REFORMER IS',F7.2/J 

109 FORMAT Cl’,' THE STEAM/METHAN'OL RATIO IN THE REFORMER IS',F7.2/> 

110 FORMAT Cl’,' THE STEAM/NAPHTHA RATIO IN THE REFORMER IS',F7.2/) 

111 FORMAT (IX, 'THE DUTY OF HEAT EXCHANGER (CAL.) (0 MEANS NO THIS NO.- 
1 HEAT EXCHANGER)' ) 

112 FORMAT (IX, 'THE SURFACE AREA OF HEAT EXCHANGER (M**2) (0 MEANS NO - 
1 THIS NO. HEAT EXCHANGER)') 

113 FORMAT (IX, 'THE EFFICIENCY OF HEAT EXCHANGER (0 MEANS NO THIS NO. - 
1 HEAT EXCHANGER OR IS CONDENSER)') 

114 FORMAT (IX, 'THE POWER OF AIR COMPRESSOR :' ,F10 . 2 ,' HP '/IX, 'THE PCWER- 
1 OF METHANE COMPRESSOR :' ,F10 . 2 , 'HP V1X, ' THE POWER OF PUMP :',F10.5- 
2, 'HP' ) 

115 FORMAT (IX, 'THE POWER OF AIR COMPRESSOR :', F10 . 2 ,' HP V1X ,' THE POWER- 
1 OF METHANOL COMPRESSOR :', F10 . 5 , ‘HP V1X ,' THE POWER OF PUMP :',F10.- 
25, 'HP' ) 

116 FORMAT (IX, 'THE POWER OF AIR COMPRESSOR :', F10 . 2 , ‘ HP V1X, ' THE POWER - 
1 OF NAPHTHA PUMP : ' , F10 . 5 , ' HP V1X , ' THE FOWER OF PUMP : ' , F10 . 5 , ' HP 1 ) 

117 FORMAT ( //IX, ' THE AC OUTPUT IS ' , F7. 2 , ' KW '// ) 

118 FORMAT! ' THE FUEL CELL IS OPERATING UNDER THESE CONDITIONS ',6X,/'- 
1 THE OPERATING TEMPERATURE :',F7.2,'KV THE OPERATING FRESSURE: ',- 
2F5 . 2 , ' ATM ' / ' THE OPEN CIRCUIT POTENTIAL: ', F8 . 3 , ' VV 

3' THE OPERATING POTENTIAL: ', F8 . 3 , ' VV THE CURRENT DENSITY :■ ,F8. 3- 
4,'A/CM*#2V THE CATALYST LOADING: ', F8 . 3, ' PT/CM**2 V 
5 ' THE FUEL UTILIZATION: ', F5 . 3/ 

6' THE OXYGEN UTILIZATION: ', F5 . 3 ) 

119 FORMATt '1' ) 

END 

SUBROUTINE EURN( DNS , I , I J ) 

C**«*******#******#****»****XJ <»********■***»**»»**#*«*#***#«**»*»*#»***#* 

C THIS SUBROUTINE IS TO CALCULATE THE MASS BALANCE OF BURNER 
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0021100 

0021200 
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C#* a*###****#*#**#*####*#**#***#***#*******-**#**** ***#»******»*»*******» 
C ASSUMPTION: 

C (1) THE COMBUSTION IS COMPLETE 

C X : THE AMOUNT OF 02 REACTED 

C XT : THE AMOUNT OF C02 PROOUCED 

C T : THE AMOUNT OF H20 PRODUCED 

C IJ: 1 =MATHANE AS INPUT GAS 

C 2 =METHANOL AS INPUT GAS 

C 3 =NAPHTHA AS INPUT GAS 

DIMENSION DNStI) 

C CALCULATE THE NECESSARY AMOUNT OF 02 
IF( IJ .EQ. 3 ) GO TO 4 
IFdJ.EQ.2) GO TO 2 

X=1./2.*0NS(3)+1./2.*DNS(5)+2.*0NS<1) 

GO TO 3 

2 X=.5*DNS(3) + .5«DNS(5H1.5*DNS(1) 

GO TO 3 

4 X-0 . 5*DNS( 3 )+0 . 5#DNS( 5 ) + 15 . *DNS( 1 ) 

XY=0NS(3)+7.*DNS(1> 

Y=DNS(5)+8.*DNStl) 

GO TO 5 

3 CONTINUE 
XY=DNSl 3 )+DNS( 1 ) 

Y=DNS(5)+2.«DNS(1) 

C CALCULATE THE EXIT COMPOSITION 

5 DNS(1)=0. 

DNS( 3 ) = 0 . 

DNS( 5 )-0 . 

DNS( 2 )-DNS( 2 )-X 
DNSt 4 )=DNS( 4 )+XY 
DNS( 6 )=DNS( 6 )+Y 
RETURN 
END 

SUBROUTINE CDPH ( QTI , TCI ,TCO , DNSC , P , I , qTS ) 
C*********************************************************************** 
C THIS SUBROUTINE IS TO CAL. THE ENERGY ANALYSIS FOR E-7 AND E-10 
C*********************************************************************** 
C DEFINITICN: 

C QTi: TOTAL HEAT TRANSFER FROM FUEL CELL 
C TCB= BOILING FOINT OF FIXED PRESSURE 

DIMENSION GS( 7) ,HS( 7 ) ,HCAS( 7 ) ,HCBS( 7 ) ,HCCS( 7 ) 

DIMENSION DNSC(I) 

COMMON /ETHDA/ GS ,HS ,HCAS ,HCBS ,HCCS 

com: ion/tc/tc/cons/a , b 

CCMMON/U1/ U 

C ASSUME THE SATURATED .PRESSURE IS EXP(A-B/T) 

TCB=B/(A-ALC3(P) ) 

C ASSUME THE SATURATED STEAM OUTPUT 
TCO=TCB 

C ASSUME THE LATENT HEAT CAL. BY WATSON CORRELATION 
TC1=TCB/TC 

C ASSUME THE AVERAGE HEAT CAPACITY OF WATER IS 1 CAL/G-K 

QT=< ( 1. -TCI )/( 1.-0.577 38*9700. 0*DNSC< 6 > + ( TCB-TCI )*1 . »18 . 
1*DNSCC6) 
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0025930 
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0026302 

0026309 

0026506 

C026308 

0026310 

0026312 

0026319 

0026316 

0026318 

0026320 

0026322 

0026329 

0026326 

0026328 

0026330 

0026332 

0026339 

0026336 

0026338 

0026390 

0026392 

0026399 

0026396 

0026398 

0026350 

0026352 

0026900 

0026500 

0026520 

C026700 

0026800 

0026900 
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0027100 

0027120 

0027200 

0027300 

0027900 

0027500 

0027600 

0027700 

0027800 

0027900 

0028000 

0028100 

0028200 


IF(QT.GE.QTI) GO TO 1 

QTI-QTI-QT 

QTS=QT 

GO TO 2 

1 TCO-TCB 

2 CONTINUE 
RETURN 
END 

SUEROUTINE COMPC DNS , TIN. TOUT, PIN, POUT, ROW ,GAG, I , IP ) 
C*********************************************************************** 
C THIS SUBROUTINE CALCULATES THE BALANCE OF COMPRESSOR 

C *********************************************************************** 
C ASSUMPTION: 

C (1). IDEAL GAS BEHAVIOR 

C GAG= RATIO OF HEAT CAPACITY 

C WS! SHIFT WORK 

C POW: POWER REQUIRED! HP 

C Vo: SPECIFIC VOLUME OF GAS AT APPLIED CONDITION! M#*3/G-M0LE 

DIMENSION DNS(I) 

TDNS=0. 

DO 1 IA=1.I 

IF ( QNS( IA ) . EQ . 0 . ) GO TO 1 
TDNS=TDNS+DNSCIA) 

1 CONTINUE 

IF( IP. EQ. 2 ) GO TO 2 

TOUT=TIN*( POUT/PIN )**( ( GAG-1 . J/GAG ) 

WS-GAG*1.937*TIN»1.8*( (POUT/PIN J**( (GAG-1. J/GAG )-l. J/( GAG-1. ) 
FOW=US*TDNS/691900. 

RETURN 

2 TOUT-TIN 

WS=1.987*TIN*1.8*ALOG( POUT/PIN) 

FCW=WS*TDNS/69190(J. 

RETURN 

END 

SUBROUTINE COND( THI .THO.DNSH ,QT,CH ,1 J 
C*********************************************************************** 
C THIS SUEROUTINE ESTIMATES THE HEAT DUTY IN THE CONDENSER 
C*** ******************************************************* ************* 
C DEFINITION IS THE SAME AS HEXC 

DIMENSION GS( 7 ) ,HS( 7 ) >HCAS( 7 ) ,HCBS( 7 ) >HCCSt 7 ) 

DIMENSION DNSH(I) 

COMMON /ETHDA/ GS .HS >HCAS ,HCBS . HCCS 
COMMON/TC/TC 

C CAL. THE MEAN TEMP. OF HOT SIDE 
THM=( THI+THO J/2 . 

C CAL. THE CAPACITY RATE OF FLUID OF HOT SIDE 
CH=0 . 

DO 1 IA=1,I 

IF( DNSHf IA ) . EQ. 0 . ) GO TO 1 

CH=CH+DNSH( IA )*( HCAS( IA ) +HCBS< IA )*THM+HCCS( IA )*THM**2 ) 

1 CONTINUE 
QT=0. 

DO 2 IA=1 > I 

QT=QT+DNSH( IA )*( HCAS( IA )*< THI-THO)+HCBS( IA )*( THI*THI-THO*THO ) 
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1+HCCS(IA)*(THI#*3-TH0**3) ) 

2 CONTINUE 

QTL=l (l.-( THO/TC))/( 1.-0.577 ) )**0. 38*9700. *DNSH( 6) 

QT=QT +QTL 

RETURN 

END 

SUBROUTINE CONV( XV, YV,NR ,NC ) 

C##**** **#*#***»#*****##*##***»*****»##*##****#***»»«»***»*#*#*##***»*#» 
C THIS SUBROUTINE USES WE3STEIN METHOD FOR ALGEBRAIC CONVERGENCE 
C****###****##*#*#*#**»************************************************* 
C DEFINITION: 

C XV : TRIAL VALUE 

C YV: CALCULATED VALUE 

C NC= CONVERSE INDEX 

C NC=1 CONVERGE 

C NO =2 NONCONVERGE 

DIMENSION XA( 2 ) , YA( 2 ) 

IF ( ABS( ( XV-YV )/( XV+YV ) ) . LT . 0 . 001 ) GO TO 2 
IF(NC.LE.l) GO TO 1 

XT=( XA( NR )#YV-YA(NR )*XV)/(XA(NR )-XV+YV-YA( NR ) ) 

XA(NR)=XV 
YA( NR ) = YV 
XV=XT 
RETURN 

1 XA( NR )=XV 
YA(NR)=YV 
XV=YV 
NC=2 
RETURN 

2 XV = YV 
NC=1 
RETURN 
END 

SUBROUTINE DIVID( TIN .TOUT1 ,TOUT2 ,DNS,DNS1,DNS2 ,GARM,I ) 
C********#*********#****»#**#****#*##****#**#**#*#*******#**#******##*** 
C THIS SUBROUTINE CALCULATES THE BALANCE AROUND THE DIVIDER 
C«#*****^*** a#*#*#*#***## a#*##**##****#*****##**#*#******#*##*#***##**** 
C GARM: DIVIDER FACTOR OF STREAM 32 

DIMENSION DNS( I ) ,DNS1( I ) ,DNS2( I ) 

C CAL. THE OUTLET CONDITION 
DO 1 IA=1,I 
DNS1( IA )=GARM*0NS( IA) 

1 CONTINUE 

DO 2 IA=1,I 

DNS2( IA ) = ( 1 . -GARM )#DNS( IA ) 

2 CONTINUE 
T0UT1=TIN 
TOUT2=TIN 
RETURN 
ENT) 

SUBROUTINE DMIX( DNS1 , DNS2 , DNS , TIN1 , TIN2 , TOUT , I , PIN1 , PIN2 , POUT ) 
C«****#*4#*»#*#*#**»**#*##*#*K#******#**«****## ************************* 
C THIS SUBROUTINE CALCULATES THE BALANCE AROUND THE MIXER 
C*#****#*********»*###»******«****»*«************»*#*****#**#****#*»*#** 



102 


0033600 

0033700 
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0033900 

0034000 
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0034300 
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0034500 

0034600 

0034700 

0034800 
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0036500 

0036600 

0036700 

0036800 

0036900 

0037000 

0037100 
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0038200 
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0038600 

0038700 

0038800 

0038900 


DIMENSION GS( 7) ,HS< 7) ,HCAS( 7 ) ,HCBS( 7) ,HCCS( 7) 

DIMENSION DNS1( I ) ,DNS2( I ) ,DNS( I ) 

COMMON /ETHDA/ GS . HS , HCAS , HCBS .HCCS 
C CAL. THE TOTAL THERMAL CONST. 

TCAS1=0. 

TCBS1=0. 

TCCS1=0. 

TCAS2=0. 

TCBS2=0. 

TCCS2=0 . 

DO 1 IA=1 i I 

TCAS1=TCAS1+DNS1( IA )*HCASC IA) 

TCAS2=TCAS2+DNS2( IA )*HCAS( IA ) 

TCBS1=TCBS1+DNS1< IA )*HCBS( I A ) 

TCBS2-TCBS2+DNS2( IA)*HCBS< IA) 

TCCS1=TCCS1+DNS1< IA)»HCCS( IA ) 

TCCS2=TCCS2+DNS2( IA)*HCCS( IA) 

1 CONTINUE 

C ASSUME THE INITIAL VALUE 
TOUT=( TIN1+TIN2 )/2 . 

C CAL. THE ENERGY BALANCE 

2 TOUTC=tTCASl*TINl+TCA52#(TIN2-TOUT)+TCBSl/2.*(TINl**2-TOUT*«2) 
1+TCBS2/2 . *( TIN2**2-TOUT**2 )+TCCSl*( TIN1**3-T0UT#*3 )/3. +TCCS2 
2#(TIN2**3-TOUT**3)/3. 1/TCAS1 

CALL CONV( TOUT > TOUTC > 1 , NC J 
GO TO ( 3 1 2 ) >NC 

C CAL. AND WRITE THE OUTLET COMPOSITION 

3 DO 4 IA=1,I 

DNS( IA ) =DNS1 ( IA )+DNS2( IA ) 

4 CONTINUE 
TDNS1=0. 

TDNS2=0. 

DO 5 IA=1 , 7 
TDNS1=TDNS1+DNS1(IA) 

TDNS2=TDNS2+DNS2( IA ) 

5 CONTINUE 

C ASSUME PRESSURE DROP TO BE 3X AT MIXER 

PCUT=( TDNS1+TDNS2 )/( TDNS1*TIN1/PIN1+TDNS2*TIN2/PIN2 )«TOUT*0 . 97 

RETURN 

END 

SUBROUTINE ENFUf DNSA.DNSC ,TINC ,TINA .TOUT, TOP ,POP, VOP.UT ,M, I ,QT, WE , 
lOUjCD »PT, I J ) 

C THIS SUBROUTINE IS TO CAL. ENERGY BALANCE OF FUEL CELL 

C#******#*##**#*#*»#*»*#*##**###*##**#*#*****#*»*#***#«#**#»*##**#**»«** 
C DEFINITION: 

C DGR : FREE ENERGY CHANGE AT FUEL CELL CONDITION 

C E: FUEL CELL EQU. POTENTIAL 

C EC: HEATING VALUE EFFICIENCY 

C EFC: FUEL CELL EFFICIENCY 

C El: CURRENT EFFICIENCY 

C EO: FUEL CELL STANDARD EQU. POTENTIAL 

C EV: VOLTAGE EFFICIENCY 

C FCH4: FRACTION OF METHANE 
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0039000 C FME • FRACTION OF METHANOL 

0039100 C FNAP: FRACTION OF NAPHTHA 

0039200 C M= INDEX OF OUTLET CONDITION 

0039300 C M=1 OUTLET TEMP. IS FIXED TO TOPFC 

0039900 C M=2 OUTLET TEMP. IS NOT FIXED 

0039500 C PCO: MOL. FRACTION OF CO 

0039600 C Q= TOTAL HEAT RELEASE PER HR. 

0039700 C SIS: INTEGRATION CONST. IN CALCULATING S 

C039e00 C TAUHO: MOL. FRACTION OF AVAILABLE H20 

0039900 C UT: FUEL UTILIZATION 

0090000 C VOP: ACTUAL FUEL CELL POTENTIAL 

0090100 C L'E: ELECTRICAL WORK 

0090200 C HLHVt I ) : LOWER HEATING VALUE OF GAS I 

0090300 DIMENSION GSC 7 ) ,HS( 7 ) ,HCAS( 7 ) ,HCBS( 7 ) ,HCCS( 7 ) 

0090900 DIMENSION HLHVt 7 ) ,DNSA( 7 ) ,DNSC( 7 ) ,DNS( 7 ) ,DNSCI< 7 ) ,DNSAO< 7 ) , 

0090500 1DNSAK7) 

0090600 COMMON /ETHDA/ GS,HS,HCAS,HCBS,HCCS 

0090700 COMMON/CONFC/ E , ETH , El . EV, EC.EFC 

0090900 COMMON/HLHV1/ HLHV 

0091000 C CAL. THE MOL. FRACTION OF AVAILABLE HYDROGEN 
0091100 TDNS=0. 

0091200 TDNSC=0. 

0091300 DO 1 IA=1 i I 

0091900 IFtDNSAtlAJ.EQ.O. ) GO TO 1 

0091500 TDNS=TDNS+ONSA(IA) 

C091600 IFtDNSCdAJ.EQ.O. ) GO TO 1 

0091700 TDNSC=TDNSC+DNSC(IA) 

009ie00 1 CONTINUE 

0091900 IF(IJ.EQ.l) FCH9=DNSA< 1 J/TDNS 

0092000 IF! IJ . EQ . 2 ) FME=DNSAt 1 J/TDNS 

0092100 IF! IJ . EQ. 3) FNAP-DNSA! 1 J/TDNS 

0092200 AHLU=DHSA(5)/T0NS 

0092300 TAUHO=DNSA( 6 J/TDNS 

0092900 PCO=DNSA( 3 J/TDNS 

0092500 PPH2=SQRT( DNSA! 5 J/TDNS*DNSA< 5 )*l 1-UT )/( TDNS-DNSAt 5 J#UT ) J 

0092600 PPCO=SQRT( DNSA! 3 )/TDNS*DNSA( 3 )/( TDNS-DNSAt 5 )*UT ) ) 

00927C0 PF02=SQRTt ONSCt 2 J/TDNSC*DNSC( 2 )*( 1-OU )/( TDNSC+DNSCt 2 J*OU ) J 

0092800 PFH20=SCRT( DNSCt 6 )/TDNSC*DNSC( 6 )/( TDNSC+DNSCt 2 )*OU J J 

0092900 CALL VINEWt 1 , VOP, CD ,TOP , POP, PFH2 , PP02 , PPH20,PPC0, 

0093000 1X0) 

0093100 C CAL. THE OPEN-CIRCUIT POTENTIAL 
0093200 DHCAS=HCASt 6 J-l ./2 . *HCAS( 2 )-HCAS< 5 ) 

0093300 DHCBS=HCBS( 6 J-l ./2 . *HC6S( 2 J-HCES! 5 ) 

0093900 DHCCS=HCCS( 6 J-l ./2 . *HCCS( 2 )-HCCS( 5 ) 

009350 0 DH0=HS(6 J-l ./2 .*HSt 2 J-HS! 5 )-DHCAS*298. -1 ,/2 .*DHCBS»298 .**2-l./3.* - 

0C93600 1DHCCS^29S.**3 

0093700 SIS=DHO/298. +DHCAS-( GS( 6 J-l ./2 . *SS( 2 J-GS! 5 ) J/298. -DHCAS*ALOG( 298 . )- 

,0093800 l-DHCBS*298./2.-DHCCS*298.**2/6. 

’0093900 DG=DH0+(CNCAS-5IS)*T0P-DHCAS*AL0G(T0P)*T0P-DHCBS/2.*T0P**2 

0099000 l-DHCCS/6 . *T0P**3 

00991C0 IFITOP.EQ.963. ) DG=-52798. 

0099200 EO=-DG/2./23060. 

0099300 E=EO+1.987#TOP/2./23060.*ALOG( AHLU*0 . 21**0 . 5/TAUHO*POP»*0 . 5 ) 

0099900 C CAL. FREE ENERGY CHANGE AT FUEL CELL CONDITION 
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DGR=-2. *23060. *E*AHLU*TDNS 
CAL. THE EFFICIENCY 
EV=VOP/E 
EI-UT 

CAL. THE ELECTRICAL WORK 
U'E = -EV*EI*DCR 

CAL. THE THERMODYNAMIC EFFICIENCY 

DH = ( -57798. +0HCAS*( TOP-298. )+0HCBS/2.*(T0P**2-298.**2 >+DHCCS# 
l(TOP**3-298.**3)/3. )*AHLU 
IF ( TOP . EQ . 463 ) DH = -53186 . *AHLU 
ETH=DGR/DH/TDNS 
IF (IJ.EQ.l) 

1DHC=FCH4*(HLHV< 1 ] + ( HCASt 4) + 2 . *HCAS( 6 )-HCAS( 1 )-2 .*HCAS( 2 ) ) 

1*( TOP-298. )+(HCBS(4)+2.*HCBS(6)-KCBS( 1 )-2.*HCBS(2) )/Z. 
2#(T0P**2-298.**2 >*( HCCSl 4 )+2 .»HCCS( 6 J-HCCSl 1 )-2 . *HCCS( 2 ) )/3. 

3*( TOP** 3-2 98.*# 3 ) ) 

IF (IJ.Eq.2) 

10HC=FME#( HLHV( 1 ) + ( HCAS( 4 ) + 2 . *HCAS(6 )-HCAS( 1 )-l .5*HCASl 2 ) ) 

1*( TCP-298. ) + (HCBS(4) + 2. *HCB5( 6 )-HCBS( 1 )-l .5*HCBS( 2 ) )/2. 
2#(TOP**2-298.**2)+(HCCS(4)+2.*HCCS(6 )-HCCS( 1 )-l .5*HCCS( 2 ) )/3. 

3*( TOP**3-298. **3 ) ) 

IF (IJ.EQ.3) 

10HC=FNAP*(HLHV(1)+(7.»HCAS(4J+8.*HCAS(6 )-HCAS( 1 )-15. *HCAS( 2 ) ) 

2*( TOP-293 . ) 

1+C6.*HCBSC6)+7.*HCBS(4)-HCBS(1)-15.*HCBS(2))* < TOP**2-29e.**2 )/2. 
3+ ( 8 . *KCCS( 6 )+7. *HCCS( 4 ) -HCCS( 1 ) -15 . *HCCS( 2 ) )/3 . 
3*(TOP**3-298.**3)) 

OHC=DHC+FCO*( HLHVt 3 )+( HCASt 4 )-HCAS( 3 )-0 . 5*HCAS( 2 ) ) 

4*( TOP-298 . )+( HCBS( 4 )-HCBS( 3J-1./2 .*HCBS( 2 ) )/2 . 

5#(TOP**2-298.**2 M-t HCCSt 4 )-HCCS( 3 )-l ./2 . *HCCS( 2 ) )/3.*(T0P##3-293. 
6**3 ) )+AHLU*( HLHVl 5 ) + ( HCAS( 6 ) -HCASt 5 )-l ./2 . *HCAS( 2 ) )*( TOP-298 . ) 
7+(KCES(6)-HCBS(5)-l./2.*HCBS(2 ) )/2 . *t TOP**2-298 . «*2 ) + ( HCCSt 6 ) 
8-HCCS(5)-l./2.*HCCS(2) )/3 . *< TOP**3-298 . **3 ) ) 

EC=DH/DHC 
EFC=ETH*EV*EI*EC 
DO 2 IA=1,I 
DNSCI(IA)=DNSC(IA) 

2 DNSAK IA )=0NSA( I A ) 

CALL FUCE(DNSAI,TOP,POP>DNSA,DNSC,DSO,DSN,DSHO,UT,I,PINF,PINA,IJ) 
DHIN=0. 

DO 3 IA=1,I 

DHIN=DHIN+DNSAI< IA )*( HCAS( IA )*( 298 . -TINA )+HCBS( IA )/2 .*( 298.**2 
1-TINA*#2 )+HCCS(IA)/3.*( 298 . **3-TINA**3 )) +DNSCK IA )*( HCAS( IA )* 
2(298.-TINC)+HCBS(IA)/2.*(298.**2-TINC**2)+HCCS(IA)/3.* 
3(298.*#3-TINC*#3)) 

3 CONTINUE 

DH=-57973. *( DNSAK 5 )-DNSA( 5 ) ) 

CAL. THE OUTLET COMPOSITION 
DO 4 IA=1 i I 

4 DNS( IA )=DNSA( IA )+DNSC( I A ) 

TCAS-0 . 

TCBS=0. 

TCCS-O. 

DO 5 IA=1,I 



0099900 

0050000 

0050100 

0050200 

0050300 

0050900 

0050500 

0050600 

0050700 

oosoeoo 

0050900 

0051000 

0C51100 

0051200 

0051300 

0051900 

0051500 

0051600 

0051700 

0051800 

0051900 

0052000 

0352100 

0052200 

0052300 

0052900 

0052500 

0052600 

0052700 

OC52SOO 

0052900 

0053000 

0053100 

0053200 

0053300 

0053900 

0053500 

0C53600 

0053700 

00538C0 

0053900 

0059000 

0059100 

0059200 

C059300 

0059900 

0059500 

0059600 

’0059700 

0059300 

C059900 

0055000 

0055100 

0055200 


TCAS=TCAS+DNS( IA )*HCAS( IA ) 

TCBS=TCBS+DNS( IA )*HCBS( IA ) 

TCCS=TCCS+DNS( IA)*HCCS( IA ) 

5 CONTINUE 
TOUT =TOP 

IF(M.EQ.l) GO TO 7 
C CAL. THE CUTLET TEMP. 

6 TOUTC=( -DHIN-DH-WE -TCBS/2 .«( (TOUT )**2-( 298 . )**2 J-TCCS/3. 

1*( <TOUT)#*3-(293. )**3 ) J/TCAS+298. 

CALL CCHV t TOUT , TOUTC , 1 , NC ) 

C-0 TO ( 8 > 6 ) i NC 

C CAL. THE TOTAL HEAT REJECTED PER HR. 

7 QT=-DHIN-DH-kE -TCAS*C TOUT-298 . )-TCBS/2 .*( T0UT»*2-29S . **2 ) 
l-TCCS/3 . *( TOUT**3-298 . **3 ) 

8 CONTINUE 

C 1 KU’HR=86C076CAL 
HE-KE/860076. 

RETURN 

END 

SUBROUTINE ENRE( DNS > TOP i POP) TIN i TOUT i I i IP , IJ ) 

C*#***»*S* **»#**#****#*#**)<#*»#*»»*»*»********»»#**#***»**»*****#***»*** 

C THIS SUBROUTINE CALCULATES THE ENERGY BALANCE OF REFORMER 
C##*********#**#**##******#**-#)****#***###*#****##****#*###*****#*****##* 
DIMENSION GS( 7 ) ,HS( 7 ) ,HCAS( 7 ) ,HCBS( 7 ) ,HCCS( 7 ) 

DIMENSION DNS(7),DINS(7),X(2) 

CCMHCN/ETHDA/ GS , HS ,HCAS ,HCBS ,HCCS 
C STORE THE INLET COMPOSITION 
DO 1 IA=1.I 
DINS(IA)=CNS(IA) 

1 CONTINUE 

C CALCULATE THE OUTLET COMP. 

CALL REF ( DNS) TOP) POP > X > I > I J ) 

IF ( IP . EQ . 1 ) GO TO 2 

TOUT=TIN 

GO TO 6 

2 CONTINUE 

C CALCULATE THE ENTHALPY CHANGE WITH TEMP. OF INLET GAS 
DHIN-0. 

DO 3 IA=1,I 

IF( 0INS( IA ) . EQ. 0 . ) GO TO 3 

CHIN=DHIN+DIN5( I A )*( HCAS( IA )*( 2 98 . -TIN )+HCBS( IA )/2 . *( ( 298 . )**2 
1-TIN**2 )+HCCS(IA)/3.*( (298. J**3- ( TIN )«*3 H 

3 CONTINUE 

C CALCULATE THE ENTHALPY CHANGE OF REACTION 
IF (IJ.EQ.l) 

1DH1 = (HS(3)+3.«HS(5)-HS(1)-HS(6))*X(1) 

IF (IJ.EQ.2) 

1DH1=(H3(9)+3.*HS(5I-HS( 1)-HS(6 ) )*X( 1 ) 

IF (IJ.EQ.3) 

1DH1 = (7.*HS(3) + 15.#HS(5)-HS(1)- 7 .*HS< 6 ) }*X( 1 ) 

IF( IJ . EQ. 2 ) 

1DH2 = (HS(3) + 2.*HS(5)-HS(1) )*X(2) 

IF ((IJ.EQ.l). OR. (IJ.EQ.3) ) 

10H2 = (HS(9)+HS(5)-HS(3)-HS(6) )*X(2) 
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0055300 

0055400 

0055500 

0055600 

0055700 

0055300 

0055900 

0056000 

0056100 

0056200 

0056300 

0056400 

C056500 

0056600 

0056700 

0056800 

0056900 

0057000 

0057100 

0057200 

0057300 

0057400 

0057500 

0057600 

0057700 

0057800 

0057900 

0053000 

0053100 

0056200 

C058300 

C05S400 

0058500 

005e600 

0058700 

0058800 

0058900 

0059000 

0059100 

0059200 

0059300 

0059400 

0059500 

0059600 

0059700 

0059800 

0059900 

0060000 

0060100 

0060200 

0060300 

0060400 

0060500 

0060600 


DH=0H1+DH2 

C CALCULATE THE TOTAL HEAT CAP. CONSTANT OF OUTLET GAS 
TCAS=0. 

TC6S=0. 

TCCS=0. 

DO 4 IA = 1 i I 

IF ( DNS( IA ) . EQ. 0 . ) GO TO 4 
TCAS=TCAS+DN5( IA >*HCAS< IA ) 

TCBS=TCBS+DNS( IA )*HCBS( IA ) 

TCCS=TCCS+DNS< IA )*HCCS( IA ) 

4 CONTINUE 

C CALCULATE THE MAX. TEMP. OF OUTLET GAS 
TOUT=TOP-500. 

5 T0UTC=!-0H-DHIN-TCBS/2.*( (TOUT) **2-< 298. )#*2 )-TCCS/3.*( ( TOUT )#*3 
l-( 298. )»*3 ) J/TCAS+298. 

CALL CONV ! TOUT . TOUT C . 1 i NC ) 

GO TO ( 6 > 5 ) i NC 

6 CONTINUE 
RETURN 
END 

SUBROUTINE ENSH ( ONS , TIN , TOUT , TOP , PIN , I , IP , I J ) 
C***#******#*#**####*****#*****#****«***»**#*****##*#*»*************»*»* 
C THIS SUBROUTINE CALCULATES THE ENERGY BALANCE OF SHIFT CONVERTER 
C*#*####*####*###*#*###*##»*»***#**#****##***##*###**»»#***«*#***#*#»»** 
DIMENSION GS( 7 ) ,HS( 7 ) ,HCAS( 7 ) ,HCBS< 7 ) ,HCCS( 7 ) 

DIMENSION DNS( 7 ) >DINS( 7 ) 

CCMMON/ETHD A/ GS , HS , HC AS , HCBS , HCCS 
C STORE THE INITIAL COMPOSITION 
DO 1 IA=1,I 

1 DINSt IA )=DNS( IA ) 

CALL PDSHtDINS, PIN, POUT, TOP, l.IJ) 

FOP=( PIN+POUT )/2 . 

C CALCULATE THE OUTLET COMPOSITION 
CALL SHIFT! DNS , TOP , POP , X, I ) 

IF! IP. EQ. 1 ) GO TO 2 

TOUT=TIN 

GO TO 6 

2 CONTINUE 

C CALCULATE THE ENTHALPY CHANGE WITH TEMP. OF INLET GAS 
DHIN=0. 

DO 3 IA=1,I 

IF ( ONS! IA ) . EQ . 0 . ) GO TO 3 

DHIN=DHIN+DINS(IA)#( HCAS! IA )*( 298. -TIN J+HCBS! IA )/2 .#( ( 298. )**2- 
1! TIN )»*2 ) +HCCS! IA)/3.#!(293. )«*3-< TIN )#*3 ) ) 

3 CONTINUE 

C CALCULATE THE ENTHALPY CHANGE OF REACTION 
DH=(HS(4)+HS(5)-HS( 3)-HS(6))*X 
TCAS=0. 

TCBS = 0 . 

TCCS=0. 

DO 4 IA=1,I 

IF ( DNS! IA ) . EQ. 0 . ) GO TO 4 
TCAS-TCAS+DNS! IA )#HCAS( IA ) 

TCBS=TCBS+DNS( IA )*HC8S( IA ) 
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0060700 

0060800 

0060900 

0061000 

0061100 

0061200 

0061300 

0061400 

0061500 

0061600 

0061700 

0061800 

0061900 

0062000 

0062100 

0062200 

0062300 

0062400 

0062500 

0062600 

0062700 

0062800 

0062900 

0063000 

0063100 

0063200 

0063300 

0063400 

0063500 

0063600 

0063700 

00638C0 

0063900 

0064000 

0C64100 

0064200 

0064300 

0064400 

0064500 

0064600 

0064700 

0064300 

0C64900 

0065000 

0065100 

0065200 

0065300 

.0065400 

"0065500 

0065600 

0C65700 

0065300 

0065900 

0066000 


TCCS=TCCS+DNS( 1A )*HCCS( IA ) 

4 CONTINUE 

C CALCULATE THE MAX. TEMPERATURE OF OUTLET GAS 
TOUT=TOP+10. 

5 TOUTC=( -DH-DHIN-TC8S/2 . *( ( TOUT )**2-( 298. )**2 )-TCCS/3.«( ( T0UT)**3 - 

l-( 298. )**3) J/TCAS+298. 

CALL C0NV( TOUT >TOUTC >1 >NC ) 

GO TO ( 6 , 5 ) > NC 

6 CONTINUE 
RETURN 
END 

SUBROUTINE EQUK(NNS,TOP,SK,I > 

C**K»* if****###*#####******##**##**#**#**####*********#****#*****#***##** 
C THIS SUBROUTINE CALCULATES THE EQUALIBRIUM CONSTANT 

C******************* *************** ************************************* 
C R: GAS CONSTANT JG-CAL/G-MOLE-K 
C TST: STANDARD TEMPERATURE; K 
C 

DIMENSION GS( 7 ) , HS( 7 ) ,HCAS( 7 ) , HCBS( 7 ) ,HCCS( 7 ) 

DIMENSION NNS(I) 

COMMON /ETHDA/ GS,HS.HCAS,HCBS,HCCS 
DATA R/1.987/ 

DATA TST/298./ 

C CAL. THE TOTAL HEAT CAPACITY CONSTANT 
TCAS-0 . 

TCBS=0. 

TCCS=0. 

DO 1 IA=1 > I 

IF( NNS( IA ) . EQ. 0 ) GO TO 1 
TCAS=TCAS+NNS( IA)*HCAS( IA ) 

TCBS=TCBS*NNS( IA ) *HCBSl IA ) 

TCCS=TCCS+NNS( IA )*HCCS( IA ) 

1 CONTINUE 

C CAL. HEAT CHANGE OF REACTION 
OH-O. 

DO 2 IA=1,I 

IF ( NNS( IA ) . EQ . 0 ) GO TO 2 
DH=DH+NNS(IA)*HS(IA) 

2 CONTINUE 

C CAL. FREE ENERGY OF REACTION 
DG=0. 

DO 3 IA=1 > I 

IF( NNS( IA ) . EQ. 0 ) GO TO 3 
DG=DG+NNS(IA)#GS(IA) 

3 CONTINUE 

C CAL. HEAT CONST. 

DHQ=DH-TCAS#TST-TCBS*TST**2/2 . -TCCS*TST**3/3 . 

C CAL. CONST. AI 

AI= ( DHO-CG-TCAS*TST*ALOG( TST ) -TCBS/2 . *TST*#2-TCCS/6 . *TST**3 l/TST/R 
C CAL. ECU. CONST. 

SK=EXP( -DHO/R/TOP+TCAS/R*ALOG( TOP )+TCBS/2 ,*TOP/R+TCCS/6 ,/R*TOP#*2 
1+AI) 

RETURN 

END 
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0066100 

0066200 

0066300 

0066400 

0066500 

0066600 

0066700 

0066600 

0066900 

0067000 

0067100 

0067200 

0067300 

0067400 

0067500 

0067600 

0067700 

0067600 

0067900 

0066000 

0066100 

0068200 

0066300 

0068400 

0068500 

0068600 

0068700 

0068300 

0068900 

0069000 

0069100 

0069200 

0069300 

0069400 

0069500 

0069600 

0069700 

0069800 

0069900 

0070000 

0070100 

0070200 

0070300 

0070400 

007C500 

C070600 

0070700 

0070800 

0070900 

0071000 

0071100 

0071200 

0071300 

0071400 


SUBROUTINE F LAME( DNS, TIN ,TF , I , I J ) 

C**#*#» a#*#******#****#**#*#*###*#**##*******#***#***#***#**#***#######* 
C THIS SUBROUTINE CALCULATES THE MAX. ADIABATIC FLAME TEMPERATURE 
C FOR 203X THEORETICAL AIR 

C*#*****# a**##*#**####*#*###**#******#**#*#*#**###*#*#*###*#****###***** 

C ASSUMPTION: 

C (1). THE COMBUSTION PROCESS GOES TO COMPLETION 

C (2). HEAT LOSSES ARE NEGLIGIBLE 

C (3). NEGLIGIBLE DISSOCIATION OF THE PRODUCTS OF COMBUSTION 

C (4). PRESSURE IS LOW AROUND 1ATM 

C TF: MAX. TEMPERATURE OF ADIABATIC FLAME 

DIMENSION GS( 7 ) ,HS( 7 ) ,HCAS( 7 ) ,HCBS( 7 ) ,HCCS( 7 ) 

DIMENSION DNS( 7 ) ,DINS( 7 ) 

CCMMON/ETHDA/ GS , HS , HCAS , HCBS , HCCS 
COMMON /EXT/EXT 
C STORE THE INLET FLUID 
DO 1 IA=1,I 
DINS( IA )=DNS( IA ) 

1 CONTINUE 

C CALCULATE THE EXIT FLUID 
CALL BURN! DNS, I , IJ ) 

C CALCULATE THE HEAT CHANGE OF REACTION AT 298 K 
DH = 0. 

DO 2 IA=1,I 

IF ( DNS( IA ) . EQ. 0 . ) GO TO 2 
DH=DH+DNS(IA)*HS(IA) 

2 CONTINUE 

DO 3 IA=1,I 

IF( DINS( IA ) . EQ. 0 . ) GO TO 3 
DH=DH-DINS(IA)#HS(IA) 

3 CONTINUE 

C CALCULATE THE ENTHALPY WITH TEMP. CHANGE 
DO 4 IA=1,I 

IF (DINSt IA ) . EQ . 0 . ) GO TO 4 

DH=DH+DINS( IA )*( HCAS( IA)#( 298. -TIN )+HCBS( IA)/2.*( ( 298. )**2-TIN**2)- 
l+HCCS(IA)/3.*( (298. )*»3-TIN#*3 ) ) 

4 CONTINUE 

C CALCULATE THE TOTAL HEAT CAPACITY CONSTANT 
TCAS=0. 

TCBS=0. 

TCCS=0. 

DO 5 IA=1,I 

IF( DNS( IA ) . EQ. 0 . ) GO TO 5 
TCAS=TCAS+DNS( IA )*HCAS( IA ) 

TCES=TCBS+DNS( IA )*HCB5( IA) 

TCCS=TCCS+DNS( IA )*HCCS( IA) 

5 CONTINUE 

C CALCULATE THE MAX. TEMPERATURE 
TF=TIN+500. 

6 TFC=( -DH-TCBS/2 . *( ( TF )**2-( 298. )**2 >-TCCS/3.*< (TF )**3-( 298. )#*3) ) - 
l/TCAS+298 . 

CALL CONV( TF ,TFC,1 ,NC ) 

GO TO ( 7 , 6 ) , NC 

7 CONTINUE 



0071500 

0071600 

00717C0 

0071800 

C071900 

0072000 

0072100 

0072200 

0072300 

0072400 

0072500 

0072600 

0072700 

0072800 
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0073300 
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0073500 
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C074200 

0074300 

0074400 
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0074740 
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0074765 
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0075700 
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0075900 

0076000 


RETURN 

END 

SUBROUTINE FUCE! DNS .TOP, POP.ONSA.ONSC ,OSO , DSN, DSHD, UT, I , PINF , PINA , - 
1IJ) 

C***«»* *##**#**»*#** »##**»#*»##**»*****»#»##**»****#***#******#**#*#*?■** 

C THIS SUBROUTINE CALCULATES THE MASS BALANCE OF FUEL CELL 
C****«*#***#'><K** *************************************************** *#«*** 

C X : CONSUMPTION OF H2 IN THE FUEL CELL UNDER UT UTILIZATION 
DIMENSION GS( 7) ,HS< 7 ) ,HCAS< 7 ) ,HC3S( 7 1 ,HCCS( 7 ) 

DIMENSION DNS! 7 ) iDNSAC 7 ) ,DNSC( 7 ) ,NNS< 7 ) 

COMMON /EXA/ EXA 
COMMON /HUMI/ WAT 
X=UT*DNS(5) 

C CAL. AND WRITE THE INLET COMPOSITION OF AIR 

1 DSO-! 1 . + EX A/1 00 . )#( 1 . /2 . *X-DHS( 2 ) ) 

DSN=DSC*3.76 

DSHO-! ( DSO+DSN )*28. 8 )*W AT/18 , 

C CAL. AND WRITE THE OUTLET COMPOSITION 
DO 2 IA=1,I 
DNSA! IA)=DNS(IA) 

2 CONTINUE 

DNSA! 5 ) =DNS( 5 )-X 
DO 3 IA=1,I 
CNSC! IA)=0. 

3 CONTINUE 

DH3C( 2 )=EXA/100 . *( 1 . /2 . *X-DNS( 2 ) ) 

DNSC! 6 )=DSHO+X 
DNSC! 7 )-DSN 

CALL FDFU(DNSA,DNSC 1 DNS,DSO,DSN,DSHO,POP,TOP,PINF,PINA,IJ) 

RETURN 

END 

SUBROUTINE HEPD! DNSA , DNSC ,THM ,TCt1 , HA , PINT ,PINS ,H ,DPJ , DP . NT ,TCO , TCI- 
1,IJ) 

C***»***##*****###**#****«#*****#*****#«**»*********»*****#**#****»***** 

C THIS SUBROUTINE CALCULATES PRESSURE DROP IN THE HEAT EXCHANGER 
C***»* #***#**#***»**********#*#*#*##**#****#*#***********:-**#**#**#****#* 
C DP : PRESSURE DROP ON THE SHELL SIDE 
C DPJ : PRESSURE DROP ON THE TUBE SIDE 
C REJ : REYNOLDS NUMBER OF TUBE SIDE 
C GS : SHELL SIDE MASS VELOCITY 
C FPRI : FRICTION FACTOR 

REAL IDT > IDS , DNSA! 7 ) >0NSC( 7 ) 

DIMENSION F LI ( 7 ) >C( 7 ) >CM< 7 ) > WM( 7 ) , F LJ( 7 ) iCMJt 7 ) >CJ( 7 ) 

COMMON /MM/ l,'M 

COMMON /HEPDT/ NP ,NR ,BSPAC ,ODT , PITCH ,CL , IDS ,IDT , FLOAR ,SURFC 
1,CLEN,S1TS2,DTH 

C»****#*#***#**»**##**###****#***#**»******»*» ************************** 
C HEAT EXCHANGER BASIS: 3/4 IN TUBE OD AND BUS 14 

C**»# X«##*#)t##**##**¥**#*#*»*#*#****###*****#K*#*##»**K#*»*****#»******* 
C CALCULATE NO. OF TUBES 

NT= HA /0.3043«»2/NP/CLEN/SURFC 
C CAL. NO. OF BAFFLES 
N3=CLEN/ESPAC 

C CAL. FREE AREA BETWEEN BAFFLES 



no 


0076100 
0076200 C 
0076300 
0076400 C 
0076500 
0076600 
0076700 C 
0076800 C 
C076900 
0077000 
0077100 
0077200 
0077300 
0077400 
0077500 
0077600 
0077700 
0077800 
0077900 
0078000 
C078100 
0078200 
0076300 
0078400 
0078500 
0078600 C 
0078700 
0073600 C 
0078900 
0079000 C 
0079100 
0079200 C 
0079300 
0079400 
0079500 C 
0079600 
0079700 
0079800 
0079900 
0030000 
0080100 
0080200 
0080300 
ocec40o 
0030500 
0080600 
0080700 
C08C800 
C080900 
0C31000 
0081100 
0081200 
0081300 
0081400 


FAREA=IDS/( ODT+CL )*CL*BSPAC 
CAL. CORRECTION FACTOR 
B0=N3+1. 

CAL. RATIO OF PITCH TRANSVERSE TO FLOW* TO TUBE DIA. 

XT=PITCH/ODT 
IF (IJ.GT.3) GO TO 2 
INSERT FLOW RATE OF EACH GAS 
SHELLSIDE CALCULATIONS 
FLI( 1 )=0NSA( 11/453.6 
FLI( 2 )=DNSA< 3 1/453.6 
FLI( 3 )=DNSA( 4 1/453.6 
FLI ( 4 )=DNSA( 6 )/453 . 6 
F LI ( 5 )=DN3A( 51/453.6 
F LI ( 6 )=DN5A( 7 1/453 . 6 
F LI ( 7 ) =DNSA( 21/453.6 

FI = FLI(1) + FLI(2) + FLI(3) + FLI(4) + FLI(5) + FLI(6) + FLI(7) 

DO 1 1=1,7 
CM( I ) = FLI( I )/FI 

1 CONTINUE 

AMW=CM( 1 )*WM( 1 )+CM( 2 )*WM( 3 )+CM( 3 )*WM14 ) +CM( 4 )*KMl 6 ) +CM( 5 )*WM( 5 ) 
1+CM( 6 )*NM( 7 1 +CM( 7 )*WM( 2 ) 

TF=( THM -273.16 1*1. 8+32. 

CALL CMASS( C , FLI , FI 1 
AMUI=VIS( C , TF > I J ) 

RHO=( AMN*PINS )/(0.7302*(TF+460. )) 

CAL. SHELL SIDE MASS VELOCITY ACROSS TUBES 
GS=FI*AMW/FAREA 

CAL. CONST. SBO (FOR STAGGERED TUBES) 

SBO=0 . 23+0 . ll/( XT-1 . 1**1. 08 
CAL. FRICTION FACTOR 

FFRI=SSO*(ODT*GS/AMUI)**(-O.15) 

CAL. PRESSURE DROP OF SHELL SIDE FLOW 

DP= BO*2 . *FPRI*NR*GS**2/32 . 1 74/3600. **2/RHO/2 116 . 2 

2 CONTINUE 

TUBESIDE CALCULATIONS 

FLJ(1)=DNSC( 11/453.6 
FLJ( 2 )=DNSC( 3 )/453. 6 
F LJ( 3 )=DNSC( 4 )/453. 6 
FLJ( 4 )=DNSC( 6 )/453 . 6 
FLJ( 5 )=0NSC( 51/453.6 
F L J( 6 )=DNSC( 71/453 . 6 
FLJ(7)=DNSC(2 1/453 . 6 

FJ = FLJ( ll + FLJI 2 l + FLJI 3 )+FLJ( 4 )+FLJ( 5 )+FLJ( 6 )+FLJ( 7 1 

DO 3 JJ=1 , 7 

CMJ( JJ)=FLJ( JJl/FJ 

3 CONTINUE 

Al.'M=CMJ(l )*WM< 1 )+CMJ< 2)*WM(31+CMJ(3)*Wm4)+CMJ(4)*Wm6) 

1+CMJ( 5 ) * S-.'M ( 5 )+CMJ( 6 )*WM( 7 )+CMJ( 7 !*WM( 2 1 
DT=DTH*( THM-TCM 1 
TW=DT +TCM 

TWF=(TW-273. 161*1. 8+32. 

TCF=( TCM-273. 161*1.8+32. 

CALL CMASSt CJ , FLJ , FJ 1 
AMUJ=VIS( C J , TCF , I J 1 



0081500 
0081600 
0081700 
0031800 
0031900 
C082000 
0082100 
0082200 
00-32300 
0082900 
0082500 
0082600 
0082700 
0082800 
0082900 
0083000 
0083100 
0033200 
0033300 
0083900 
0083500 
0083600 
0083700 
0083800 
0083900 
0069000 
0C89100 
0089200 
0039300 
0089900 
0089500 
0089600 
0089700 
0069800 
0089900 
00C5000 
OOS5100 
0085200 
0065300 
0085900 
0085500 
0085600 
00S57C0 
0085800 
0085900 
0086000 
0086100 
,0086200 
'0086300 
0086900 
0036500 
0086600 
0086 700 
0086S00 


AMUW=VIStCJ,TWF,IJ) 

IFUAMUJ.LE.O. ).OR.(AMUW.LE.O. )) GO TO 10 
TKC=THC(CJ,TCF,IJ) 

CP=HTCP(CJ >TCF ) IJ ) 

AROH-( AWM*PINT )/( 0 . 7302*( TCF+960 . ) ) 

GJ r F J*AUM /FLOAR/NT 
REJ=GJ*IDT/AMUJ 
IFIREJ.LE.2100. )GO TO 9 
SF = 0 . 096/RE J**0 . 2 
GO TO 5 
9 SF=16./REJ 

5 CONTINUE 

IF(S1TS2.GT. 0.715) GO TO 6 
CKC^O . 9*( 1 . 25-S2TS1 ) 

GO TO 7 

6 CKC=0 . 75* ( 1 . -S2TS1 ) 

7 CONTINUE 

CK1 = ( 1 . -S2TS1 )**2+CKC+0 . 5*( NP-1 . )/NP 
IF ( RE J .GT . 2100 . ) GO TO 8 
PHI=1 . 1*( AMUJ/AMUW )**0 . 25 
BI=l.+CKl*IDT*FHI/9./SF/CLEN 
GO TO 9 

8 PHI=1.02*< AMUJ/AMUW 1**0. 19 

BI=1 . +0 . 51*CK1*NP*( TW -TCM )*< AMUJ/AMUW )**0 . 28/( TCO-TCI ) 
1/(CP*AMUJ/TKC/AUM 1**0.6667 
C CAL. TUBESIDE DP 

9 DPJ= BI*2 . *SF*GJ**2*CLEN*NP/32 .17/3600 . **2/AROH/IDT/PHI/2116 . 2 
GO TO 12 

10 DPJ=0. 

WRITE (6ill) 

11 FORMATCTHE TRYING TEMP. IS BELOW THE LIMIT OF CAL. VISCOSITY’) 

12 CONTINUE 
RETURN 
END 

SUBROUTINE HEXCf THI ,DNSH ,ONSC .TCI ,THO , TCO ,qT ,MOD ,HA ,N , M , I , PT , PS ,NT 

l.IJ) 

C*********************************************************************** 
C THIS SUSROUTINE IS TO CAL. THE ENERGY ANALYSIS FOR HAET EXCHANGER 
C*********************************************************************** 
C cc: CAPACITY RATE CF FLUID ON COLO SIDE ,ONSC*CPC 

C CH= CAPACITY RATE OF FLUID OH HOT SIDE, DNSH*CPH 

C CMAX • MAX. CAPACITY RATE 

C CMIN: MIN. CAPACITY RATE 

C CPC : SPECIFIC HEAT OF COLD SIDE FLUID 
C CPH • SPECIFIC HEAT OF HOT SIDE FLUID 
C HE: HEAT EXCHANGER EFFECTIVENESS 

C QT: TOTAL HEAT TRANSFER RATE ACROSS HEAT EXCHANGER 
C GMAX: THE MAX. HEAT TRANSFER RATE ACROSS HEAT EXCHANGER 
C TCI : COLD SIDE INLET TEMPERATURE 
C TCO : COLD SIDE OUTLET TEMPERATURE 
C THI: HOT SIDE INLET TEMPERATURE 
C THO: HOT SIDE OUTLET TEMPERATURE 

C UA: OVERALL HEAT TRANSFER COEFFICIENT OF EXCHANGER 
C MOD: TYPE OF HEAT EXCHANGER 
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0006900 C 

0007000 C 

0007100 C 

0007200 C 

0007300 C 

0037400 C 

0007500 C 

0007600 C 

0007700 C 

0087000 

0087900 

0088000 

0083100 

ooeaooo 

ooe83oo 

0088400 C 

0068500 

0038600 

0068700 

0006300 

0086900 

0069000 

0089100 C 

0089200 

0089300 

0089400 

0089500 

0089600 

0089700 

0089800 

0039900 

0090000 C 

0090100 

0090200 

0090300 

0090400 

0090500 

0090600 

0090700 

0090000 

0090900 

0091000 C 

0091100 

0091200 

0091300 

0091400 

0091500 

0091600 

0091700 

0091600 

0091900 

0092000 

0092100 

0092200 


MOD = l COUNTERFLOW 
MOD = 2 CROSS FLOW 
M0D = 3 CONDENSER 
THM: MEAN TEt1P. OF HOT SIDE 
TCM = MEAN TEMP. OF COLD SIDE 

N= THE NUMBER OF HEAT EXCHANGER 
M-. THE INDEX OF INITIAL CONDITION 
M=1 TEMP. OF BOTH SIDES ARE KNOWN 

M=2 TEMP. OF HOT SIDE INLET AND COLD SIDE OUTLET ARE KNOWN 
DIMENSION GS< 7 ) ,HS( 7 ) ,HCAS( 7 ) ,HCBS( 7 ) ,HCCSC 7 ) 

DIMENSION DNSH(7),0NSC<7) 

COMMON /ETHDA/ GS , HS , HCAS , HCBS ,HCCS 
COMMON/U1/ U 
COMMON/CN1/ CN 
COMMON/HE/ HE 

ASSUME THE MEAN TEMP. AT COLD AND HOT SIDE 
IF( M . EQ . 2 ) GO TO 1 
THM=( THI+10.+TCIJ/2. 

TCM = ( THI-10 . +TCI )/2 . 

GO TO 2 

1 THM=( THI+TCO-50 . )/2. 

TCM=( TCO*2 . -200 . )/2. 

CAL. CC AND CH 

2 CC=0. 

CH = 0. 

DO 4 IA=1 > I 

IF< DNSC( IA ) . EQ. 0 . ) GO TO 3 

CC=CC+DNSC( IA )#( HCAS< I A )+HCBS( IA )*TCM+HCCS( IA >*TCM**2 ) 

3 IF ( DNSH( I A ) . EQ . 0 . ) GO TO 4 

CH=CH+DNSH( IA )*( HCASC IA )+HCBS( IA )*THM+HCCS< IA )*THM**2 ) 

4 CONTINUE 

CHOOSE THE CMAX. ,CMIN. 

IF(CC.GT.CH) GO TO 5 

CMAX=CH 

CHIN=CC 

GO TO 6 

5 CMAX=CC 
CMIN r CH 

6 CONTINUE 
HA=CNKCMIN/U 
UA=HA*U 

CAL. THE HEAT EXCHANGER EFFECTIVENESS 
IF(MOD.GE.2) GO TO 8 
IF ( ( CMIN/CMAX ) .GT .0.01) GO TO 7 
HE=1 . -EXP( -UA/CMIN ) 

GO TO 12 

7 HE=( 1 . -EXP( ( -UA/CMIN )*< 1. -CMIN/CMAX ) ) )/( 1 .-( CMIN/CMAX )«EXP( ( -UA 
1/CMIN )#( 1 . -CMIN/CMAX I ) ) 

GO TO 12 

8 IF( MOD . GT. 2 ) GO TO 11 

IF( ABS( CMIN/CMAX-1 . J.GE.0.01) GO TO 9 
HE = ( UA/CMIN )/(UA/CMIN+l. ) 

GO TO 12 

9 IFU CMIN/CMAX ).GT. 0.01) GO TO 10 



0092300 

0092900 

0092500 

0092600 

0092700 

C092300 

C092900 

0093000 

0093100 

0093200 

0093300 

0093900 

0093500 

0093600 

C093700 

0093800 

0093900 

0099000 

0099100 

0099200 

0099300 

0099900 

C099500 

0099600 

0099700 

0099300 

0099900 

C095000 

0095100 

0095200 

0095300 

0095900 

0095500 

0095600 

0095700 

0095710 

0095720 

0095730 

0095300 

0095900 

0096000 

0096100 

0096200 

0096300 

0096900 

0096500 

C096600 

,0096700 

’0096300 

0096900 

0097000 

0097100 

0097200 

0097300 


HE = 1 . -EXP( -UA/CMIN ) 

GO TO 12 

10 IF(CMAX.EG).CH) HE = 1 . -EXP( ( -CMAX/CHIN )*< 1 . -EXP( -UA/CMAX )) ) 

IF l CMIN.EQ.CH ) HE=( CMAX/CMIH )*< 1 . -EXP( ( -CMIN/CMAX )*( 1 . -EXPl -UA/ 
1CMIN)) ) ) 

GO TO 12 

11 HE=1.-EXP( -UA/CMIN) 

C CAL. THE OUTLET CONDITION AND TOTAL HEAT TRANSFER RATE 

12 IF( M . EQ . 2 ) GO TO 13 
THO-THI-HE*C CNIN/CH )*ITHI-TCI ) 

TCO=HE*( CMIN/CC )*( THI-TCI )+TCI 
QT=HE*CMIN*< THI-TCI) 

GO TO 19 

13 TCI= ( HE*( CMIN/CC )*THI-TCO )/( HE#( CMIN/CC )-l . ) 

THO=THI-HE*< CMIN/CH )*( THI-TCI ) 

QT=HE*CMIN*( THI-TCI) 

19 IF ( ( A3S( ( Th'O+THI )/2 . -THM ) . LT. ( ABS( ( THO+THI )/2 . +THM )*0 .005 ) 

1 ) . AND . ( ABS( ( TCO+TCI )/2 . -TCM ) . LT . I ABS( < TCO+TCI )/2 . +TCM )*0 . 005 ) ) ) 

2 GO TO 15 
THM=( THO+THI )/2 . 

TCM= ( TCI+TCO )/2 . 

GO TO 2 

15 THM=( THO+THI )/2. 

TCM=< TCO+TCI )/2. 

CALL HEFD(DNSH,DNSC,THM,TCM,HA,PT,PS,N,DPJ,DP,NT,TCO,TCI,IJ) 

IFIM.EQ.2) GO TO 16 

PT^PT-DPJ 

PS=PS-DP 

GO TO 17 

16 PT=P.T+DPJ 
PS=PS-DP 

17 CONTINUE 
RETURN 
END 

SU3R0UTINE PDFU( DNSA ,ONSC ,DNS , DSO , DSN , DSHO , POP , ATMP , PINF , PINA , I J ) 
C****** ******************************************** ************* ******** 
C THIS SUBROUTINE CALCULATES PRESSURE DROP IN THE FUEL CELL STACK 
C*#****# a*#*#*##*#**##**#**********#*#***#**#********* a**##*#****#*#**** 
REAL L( 2 ) >DP(9 ) 

DIMENSION FLI( 7) >C( 7) iCM( 7 ) »kTI( 7) >DNSA( 7) , DNSC( 7) >DNS( 7 ) .DNSIl 7 ) 
DIMENSION NTA(2), WIDA(2).D(2) 

COMMCN/FDFUT/ NTA , L , WIDA , NP 
COMMON /WM/ WM 

C**»»*#**## **#****#*»*##**##***#»«#*«***#****#**#****#***#*************# 
C BASIS: NO. 522 STACK 

C»**x*************************»##»**********#*«**********»*«************ 

IT=1 

DO 1 IA=1,7 
1 DNSI(IA)=0. 

DNSIl 9 ) =DSHO 
CNSI( 6 )-0SN 
DNSK7)=DS0 

C CAL. THE PRESSURE DROP OF FUEL SIDE 
JJ=1 
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0097400 

0097500 

0097600 

0097700 

0097500 

0097900 

0093000 

00931C0 

0090200 

0093300 

0098400 

009S500 

0098600 

0098700 

0093800 

0096900 

0099000 

0099100 

0099200 

C099300 

0099400 

0099500 

0099600 

0099700 

0099300 

0099900 

0100000 

0100100 

0100200 

0100300 

0100400 

0100500 

0100600 

0100700 

0100600 

0100900 

0101000 

01011C0 

0101200 

01013C0 

0101400 

0101500 

0101600 

0101700 

0101800 

0101900 

0102000 

0102100 

0102200 

0102300 

0102400 

0102500 

C102600 

0102700 


D( 1 )=MIDA( 1 ) 

FLI( 1 )=DNS( 1 1/453.6 
FLI( 2 )=DNS( 3 1/453.6 
FLI( 3 )=DNS( 4 ) /453.6 
FLI(4) = DNS(6) /453.6 
FLI(5)=DNS(5) /453 . 6 
FLI(6)-DNS (71/453. 6 
F LI ( 7 1 =DNS( 2 1/453 . 6 
GO TO 3 

2 FLI( 1 )=CNSA( 1 1/453.6 
FLI( 2 ) = DMSA( 3 1/453.6 
F LI ( 3 )=DNSA( 41/453.6 
F LI ( 4 1 =DNSA( 6 1/453.6 
FLI( 5 )=0NSA( 5 1/453.6 
FLI( 6 )=DNSA( 7 1/453.6 
FLI(7)=DNSA(2 1/453.6 

3 CONTINUE 
TK = ATI*1P 

FI=FLI(1) + FLI(2) + FLI(3)+FLI(4)+FLI(5) + FLI(6) + FLI(7) 

DO 4 1=1,7 
Cmi) = FLI(I)/FI 

4 CONTINUE 

ANN=CN( 1 )*HM( 1 ) +CH( 2 )*WM( 3 1 +CM( 3 )*KM( 4 1 +Ct1( 4 )*W11( 6 1 +Cf1( 5 )*WM( 5 1 
1+Ct1( 6 lawn 7)+Cf1( 7 )*KM( 2 1 
G=FI*AN'.J/HIDA( IT )**2/NP/NTA( IT 1 
TF=(TK-273. 16) *1.8+32. 

CALL Ct1ASS(C,FLI,FI) 

ANUI=VIS( C iTF i IJ 1 

RHC= ( ANN* POP )/( 0 . 7302*( TF+460 . 1) 

RE=D( IT )*G/AMUI 
CONS=57. 2 
FRIC=CCNS/RE 

DP( JJ)=RHO*(G/RHO )**2/2.*( 0.5+1. +0 . 6+FRIC*L( IT l/D (IT) 1/2116.2 
1/3600. **2/32. 174 
JJ=JJ+1 

IFUJ.EQ.2) GO TO 2 
IF( JJ . EQ. 3) GO TO 5 
IF(JJ.Eq,4) GO TO 6 
IFUJ.EQ.5) GO TO 8 
CAL. THE PRESSURE DROP OF AIR SIDE 
INSERT THE FLOW RATE OF EACH GAS 

5 FLI( 1 )=DNSC( 1 1/453.6 
FLI( 2 )=CNSC( 3 1/453.6 
F LI( 3 )=DNSC( 4 1/453 . 6 
FLI( 4 )=DN5C( 6 1/453 . 6 
FLI ( 5 1 =DNSC( 5 1/453.6 
F LI( 6 )=DNSC (71/453.6 
F LI ( 7 )=DNSC( 2 1/453 . 6 
IT=IT+1 

D ( 2 1 =NIDA( 2 ) 

GO TO 3 

6 CO 7 IA=1,7 

7 FLI(IA) = DNSI(IA 1/453.6 
GO TO 3 



0102000 

0102900 

0103000 

0103100 

0103200 

0103300 

0103900 

0103500 

0103600 

0103700 

0103300 

0103900 

0109000 

0109100 

0109200 

0109300 

0109900 

0109500 

0109600 

0109700 

0109300 

0109900 

0105000 

0105100 

0105200 

01C5300 

0105900 

C1C5500 

0105600 

0105700 

0105600 

0105900 

0106000 

0106100 

0106200 

0106300 

0106900 

0106500 

0106600 

0106700 

0106800 

0106900 

0107000 

0107100 

0107200 

0107300 

0107900 

0107500 

'0107600 

0107700 

0107000 

0107900 

0103000 

0108100 


8 CONTINUE 

PINF = PINF-( DP( 1 )+DP( 2 1 1/2 . 

PINA=PINA-CDP( 31+DPC 9 ) )/2 . 

RETURN 

END 

SUBROUTINE POSH ( DNS , PIN , POUT , TK , JK , I J 1 

CX*XXXXXXXXXXXXXXXXXXXX*XXXXXXXXXX*XXXXX*X*#XXX»XXXXXXXXXX**X»*XX*X**XXX 

C THIS SUBROUTINE CALCULATES PRESSURE DROP IN THE SHIFT CONVERTER 

CXX*XXXXXXXXXX*XX*XXXXXXXXXXXXX*XXXXXXXXKXXXXXXXX*X*XXttXXXXXXXXXX«XX«X*X 

DIMENSION D( 2 ) , AHRNC 2 ) , API 2 ) ,CLEN( 2 ) ,NT( 2 ) 

DIMENSION FLI( 7 ) ,C< 7 ) ,CM( 7 ) , WMC 7 ) ,DNS< 7 1 
CCMMON/PDSHT/ D , AHRN , AP , CLEN , NT 
COMMON /WM/ MM 
C JK=1 SHIFT CONVERTER 

C JK=2 REFORMER FOR METHANOL AND NAPHTHA FUEL 
TF=( TK-273. 16 1*1 .8+32 . 

DP=6 .»( 1 . -AHRN( JK ) )/AP( JK ) 

F LI ( 1 )=DNS( 1 )/953.6/NT( JK ) 

FLIC 2 )=DNS( 3 )/953 .6/NT( JK ) 

FLIC 3 )=DNS( 9 )/953 .6/NTC JK ) 

FLIC 9 )=DNS( 6 )/953 .6/NTC JK ) 

F Lit 5 )=DNS( 51/953. 6/NTC JK) 

FLIC 6 )=DN5( 7 1/953. 6/NTC JK 1 
FLIC 7)=0NS( 2 1/953. 6/NTC JK 1 

FI = FLI(1) + FLI(2) + FLI(3) + FLI(9) + FLIC5) + FLI(6) + FLI(7) 

DO 1 1=1,7 
CMC I ) = FLI( I l/FI 
1 CONTINUE 

AMW=CM( 1 1XWMC 1 l+CMC 2 IxwriC 3 1+CMC 3 )«WM( 9 l+CMC 9 )XKM( 6 l+CMC 5 )xWM( 5 1 
1+CMC 6 )*WM( 7 1 + CMC 7 )*WM( 2 1 
G=FI*AMM*9./( 3. 19159*D( JK )«X2 ) 

CALL CMASSC C , FLI ,FI 1 

AM'JI=VIS( C >TF , IJ 1 

RHO=( AMW*PIN )/( 0 . 7302*( TF+960 . 1 1 

DELP=CLEN( JK )X( 1 . -AHRNC JK 1 l/AHRNC JK )xx3xGXXg/DP/9 . 16E + 08 
1/RHOXC 150. x( 1. -AHRNC JK) 1XAMUI/DP/G+1. 751/2116. 2 
PO'JT=PIN-DELP 
RETURN 
END 

SUBROUTINE PUMPC DNS , TIN , TOUT , PIN , POUT , POW , 1 1 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxx*xxxxxxxxxxx*xxxxxx*x*xxx**xxxxxxxxxxx»xxx* 
C THIS SUBROUTINE IS TO CAL. THE BALANCE OF PUMP FOR WATER 
Cx***xxxxxxx*xx*xxxxxx*x*xx#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*xxxx*xxx* 
C ASSUMPTION AND DEFINITION IS THE SAME AS COMPR 
DIMENSION DNSC 7) 

DIMENSION SVC 3 1 > WMC 7 1 
COMMON/WM/ WM 
COMMON/SV/ SV.SVW 
TOUT=TIN 
C CAL. THE WORK 

C ASSUME AVERAGE SPECIFIC VOLUME OF WATER IS 0.0162 FTXX3/LBM 

P0W=SVWX199.X5. 05051X0 . 0000001*WM( 6 )X 19 . 7x( POUT-PIN )xDN5( 6 1 
1/953.6 
RETURN 
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0108200 

0108300 

0108400 

0103500 

0108600 

0103700 

0108800 

0103900 

0109000 

0109100 

0109200 

0109300 

0109400 

C1095C0 

0109600 

0109700 

0109800 

0109900 

0110000 

0110100 

0110200 

0110300 

0110400 

0110500 

0110600 

0110700 

0110800 

0110900 

0111000 

0111100 

0111200 

0111300 

0111400 

0111500 

0111600 

0111700 

0111800 

0111900 

0112000 

0112100 

0112200 

0112300 

0112400 

0112500 

0112600 

0112700 

0112800 

0112900 

0113000 

0113100 

0113200 

0113300 

0113400 

0113500 


END 

SUBROUTINE PUP< DNS , TIN , TOUT , PIN , POUT , POW , I , I J ) 
C*#»####*#***#*##**********#*#*****#»***###*#*# ##***«»*##»***##*#*#*##*» 
C THIS SUBROUTINE CALCULATES POWER NEEDED IN THE FUEL PUMP 
C**#**#*#»#*# #»***»»*# ************ ****** a***#***#*##*#****##*##**#**##** 
DIMENSION DNS( 7 ) iKM( 7) iSV( 3 ) 

CCMMON/SV/ SV.SV14 
CCMMON/Wri/ WM 
T0UT=TIN 

POW=SV( I J )*144 .*5 . 05051*0 . 0000001*WT1( I J >*14. 7*( POUT-PIN ) 

1*DNS( 1 J/453.6 
RETURN 
END 

SUBROUTINE REF ( DNS , TOP , FOP , X , I , I J ) 

C**#********#*#*»*****#*»»*#*****#******#**»»*##********»***«#***»*##*** 
C THIS SUBROUTINE IS TO CALCULATE THE MASS BALANCE OF REFORMER 
C*#*#*#*x#»***#**«»*********#*************#*#**#**####***#****##***#*#*# 
DIMENSION NNS1( 7 ) ,KNS2< 7 ) ,DNSC 7 ) ,SK1 ( 2 ) ,X( 2 ) 

C CALCULATE THE EQU. CONSTANT OF REACTION 1 
DO 1 IA=1,I 

1 NNS1 ( IA ) = 0 

IF (IJ.E0.3) GO TO 2 
IF(IJ.EQ.l) NNS1(3)=1 
IF (IJ.EQ.2) NNS1 ( 4 ) =1 
NNS1( 1 )=-] 

NNS1(6 ) = -l 
NNS1 (51=3 
GO TO 3 

2 NNS1(1)=-1 
NNS1 ( 6 )=-7 
NNS1( 3) =7 
NNSK 5) =15 

3 CALL EQUK(NNS1,T0P,SK,I> 

IHUI=0 

DO 4 IA=1,I 

4 IHUI=IHUI+NNS1(IA) 

SKI ( 1 )=SK*POP**(-IHUI) 

C CALCULATE THE EQU. CONSTANT OF REACTION 2 
DO 5 IA=1,I 

5 NNS2(IA)=0 

IF (IJ.EQ.2) GO TO 6 

NNS2( 4 ) = 1 

NNS2(5)=1 

NNS2( 3 ) = - 1 

NNS2 (6 ) = — 1 

GO TO 7 

6 NNS2(3)=1 
NNS2( 5 )=2 
NNS2( 1 ) = -l 

7 CALL EQUK( NNS2 > TOP < SK > I ) 

IHUI=0 

DO 8 IA=1,I 

8 IHUI=IHUI+NNS2(IA) 

SK1(2)=SK*P0P*«(-IHUI) 



0113600 

0113700 

0113600 

0113900 

0114000 

01141CO 

0114200 

0114300 

0114400 

0114500 

0114600 

£114700 

0114600 

0114900 

0115000 

0115100 

0115200 

0115300 

0115400 

0115500 

0115600 

0115700 

0115600 

0115900 

0116000 

0116100 

0116200 

0116300 

0116400 

0116500 

0116600 

0116700 

0116600 

0116900 

0116920 

0116940 

0116960 

0116980 

0117000 

0117020 

0117040 

0117060 

0117080 

0117100 

0117120 

0117140 

0117160 

0117180 

"0117200 

0117220 

0117240 

0117260 

0117280 

0117300 


C CALCULATE THE EXIT AMOUNT OF GAS I 
C INITIAL CONDITION 

IF IIJ.EQ.2) GO TO 10 
IF (IJ.EQ.l) GO TO 9 

C BECAUSE CF COMPUTATION PROBLEM! OVERFLOW ) NAPHTHA INPUT FUEL USINS THE 
C REASONABLE ASSUMPTION OF CONVERSION 
X( 1 )-0 . 999«DN5( 1 ) 

X(2)=2.9*0NS(1) 

DNS( 1 )=DNS( 1 )-X( 1 ) 

0NS( 3 ) =DNS( 3 ) +7 . *X(1 )-X( 2 ) 

DHS( 4 )=DNS( 4 )+X( 2 ) 

DNS l 5 )=DNS( 5 )+X( 2)+15.#X(l) 

DNS(6)=DNS(6)-X(2)-7.*X<1) 

GO TO 11 

9 Xll)=0.8*DNS(l) 

X( 2 )-0 . 35*DNS( 1 ) 

CALL SNAE(X,2,DHS,SK1,I,IJ) 

ONS( 1 1=DNS( 1 )-X( 1 ) 

DNS( 3 )=DNS( 3 )+X( 1 )-X( 2 ) 

DNS( 4 ) -DN5( 4 ) +X( 2 ) 

DHS(5)=DNS(5HX(2)+3.*X(1) 

DNS( 6 )=DNS( 6 ) —X t 2 ) — XC 1 ) 

C-0 TO 11 

10 X(1) = 0.96*DNS(1) 

X(2)=0.04#DN5ll) 

CALL SNAEC X > 2 > DNS > SKI , I >IJ ) 

DNSl 1 )=DNS( 1 )-X( 1 ) — Xt 2 ) 

DNS( 3 )=DNS( 3 HX< 2 ) 

DHS(4)=DNS14)+X( 1) 

DNS( 5 )=DNS( 5 ) + 2 . *X( 2 )+3.#X( 1 ) 

DNS( 6 )=DNS( 6 I — Xt 1 ) 

11 CONTINUE 
RETURN 
END 

SUBROUTINE SE PAR ( TIN , POP , TOUTV , TCUTL , DNS , DNS L , DNSV , I ) 

C**# ********** **************** **************************** *** *********** 
C THIS SUBROUTINE IS TO CAL. THE MASS BALANCE AROUND THE LIQUID SEPARATO 
C*********************************************************************** 
C ASSUMPTION : 

C (1). ONLY WATER EXIST IN LIQUID PHASE 

C PSAT: SATURATE FRESSURE AT T; EXPtA-B/T) FOR WATER 
C DK= EQU. CONST. OF LIGEID-VAFQR 
C xw: AMOUNT OF WATER IN LIQUID PHASE 
DIMENSION DNS( I ) ,DNSV( I ) ,DNSL( I ) 

COMMON /PS/ PS , TC/CONS/ A,B 
C CAL. THE EQU. CONST. OF LIQUID-VAPOR 
P3AT=EXPC A-B/TIN ) 

DK=PSAT /POP 

C CAL. THE EQU. AMOUNT OF LIQUID-VAPOR WATER 
TDNS r 0 . 

DO 1 IA=1,I 
TDNS=TDNS+DNS(IA) 

1 CONTINUE 

XW=(TDNS*OK-DNS(6!)/(OK-1. ) 
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0117320 

0117340 

0117360 

0117380 

0117400 

0117420 

0117440 

0117460 

0117480 

0117500 

0117520 

0117540 

0117560 

0117580 

0117600 

0117620 

0117640 

0117660 

0117680 

0117700 

0117720 

C117740 

0117760 

0117730 

0117800 

0117620 

0117840 

0117860 

0117880 

0117900 

0117920 

0117940 

0117960 

0117980 

0118000 

0118020 

0118040 

0118060 

0118080 

0118100 

0118120 

0118140 

0118160 

0118180 

0118200 

0118220 

0118240 

0118260 

0118280 

0118300 

0118320 

0116340 

0118360 

0118380 


DO 2 IA-1 > I 
DNSL(IA) = 0. 

2 CONTINUE 
0NSL(6)=XW 
DO 3 IA=1,I 
DNSVf IA)=DNS(IA) 

3 CONTINUE 

DNSVf 6 )-DNS(6 )-DNSL( 6 ) 

TCUTL=TIN 

TOUTV=TIN 

RETURN 

END 

SUBROUTINE SHIFTf DNS , TOP , POP , X , I ) 

C******#*******#***#*#**#**»*»*******#*******»#***»****#»*************** 

C THIS SUBROUTINE IS CALCULATE THE MASS BALANCE OF SHIFT CONVERTER 
C#*#*#»***#****#**##***##»*****#*********«*******#****#**#*****#*****»»* 

C ASSUMPTION: 

C (1). ONLY ONE REACTION (C0+H20_H2+C02 ) DOMINATE 

DIMENSION NNSf 7 ) >0NS( 7 ) 

F( X )=( D4+X )»( D5+X )-SK*( D3-X)#( D6-X) 

DF(X) = (D4+X) + (D5+X) + SK»( ( D3-X ) + 1 D6-X ) ) 

D1=DNS(1) 

D2=DNS(2) 

D3=DNS(3) 

D4=0NS(4> 

D5=DNSf 5) 

D6=DNS(6) 

D7-DNS( 7 ) 

C CALCULATE EQUALIBRIUM CONSTANT 
DO 1 IA=1,I 

1 NNSf IA )=0 
NNSf 4 )=1 
NNSf 5 ) = 1 
NNSf 3 ) = — 1 
NNSf 6 ) = — 1 

CALL EQUKfNNS>TOPiSK,I ) 

C NEWTON-KAFHSON METHOD TO SOLVE NONLINEAR ALGEBRAIC EQUATION 
X=0.5*DNS(3) 

DO 2 IA=1,500 
DX=ABS( F(X)/DF( X ) ) 

X=X-F(X)/DF(X) 

IF(DX.LT.O.Ol) GO TO 3 

2 CONTINUE 

3 CONTINUE 

C CALCULATE THE EXIT AMOUNT OF GAS I 
DNSf 4 l = 0NSf 4 )+X 
DNSf 5 )-DNS( 5 !+X 
OKS ( 3 ) r DNS( 3 )-X 
DNSf 6 )=DNS( 6 )-X 
RETURN 
END 

SUBROUTINE SNAE ( XY , IX , DNS , SKI , I , I J ) 

C***#***###*****#*#####»*#**####*#«***#*#***#***»*#*»******************* 

C THIS SUBROUTINE IS USING NEUTON-RAFHSON ITERATION TO SOLVE TWO NONLINE . 



0118400 

0118420 

0118440 

0118460 

0116480 

0118500 

0118520 

0118540 

0118560 

0118580 

0118600 

0118620 

0118640 

0118660 

0118680 

0118700 

0113720 

0118740 

0118760 

0118780 

onesoo 

0116820 

0118840 

0118860 

0116880 

0118900 

0118920 

0118940 

01ie960 

0118980 

0119000 

0119020 

0119040 

0119060 

0119080 

0119100 

0119120 

0119140 

0119160 

0119180 

0119200 

0119220 

0119240 

0119260 

0119280 

0119300 

0119320 

0119340 

'0119360 

0119380 

0119400 

0119420 

0119440 

0119460 


C ALGEBRAIC EQUATIONS IN REFORMER 

C******************************************** *************************** 
DIMENSION XY( IX ) iDNS( I ) >SK1< IX) 

Fl(X,Y) = (D3+Y-X)*(D5+X+3.*Y)**3-SKA*l TDNS + 2. *Y)**2*( Dl-Y )*(D6-X-Y) 
F2(XiY)=(D4+X)*(05+X+3.*Y )-SKB*( D6-X-Y )*( D3-X+Y ) 

DFX1(X>Y)=-(D5+X+3.*Y)**3+3.*(D3+Y-X )*(05+X+3.*Y)**2+SKA*(TDNS+2.- 
1*Y)**2*(D1-Y) 

DFY1( X. Y )=( D5+X+3 -*Y )**3+9.*( D3+Y-X )*< D5+X+3.*Y )**2-SKA*( -( TDNS+2 
1*Y )**2*( Dl-Y )-( TDNS+2 .*Y )**2*( D6-X-Y ) +4, *( TDNS+2 . *Y)*(D1-Y)*(D6-X - 
2-Y ) ) 

DFX21 XiY)=(05+X+3.*Y ) + ( D4+X )+SKB*( ( D3-X+Y )+(D6-X-Y ) ) 
DFY2(X,Y)=3.*(04+X)-SKB*< ( 06-X-Y ) - ( 03-X+Y ) ) 

F3(X.Y)=(D4+X)*(D5+3.*X+2.*Y)**3-SKA*( TDNS+2. *X+2 . *Y )**2*( Dl-X-Y - 
1)*(D6-X) 

F4(X,Y) = (D3+Y)*l 05+3. *X+2.*Y)**2-SKB*l TDNS+2. *X+Z . *Y )*( Dl-X-Y ) 
0FX3(X,Y)=CD5+3.*X+2.*Y)**3+9.*(D5+3.*X+2.*Y)**2*(D4+X)-2.*SKA 
1*( TDNS+2. *X+2.*Y)*( Dl-X-Y ]*< D6-X )*2 . +SKA*( TDHS+2 . *X+2 . *Y ) 
2**2*(D6-X)+SKA*(TCNS+2.*X+2.*Y)**2*( Dl-X-Y) 

DF Y3( X, Y ) = ( D4+X )*6.*( D5+3 . *X+2 . *Y )**2-4 ,*SKA*( TDNS+2. *X+2.*Y)» 

1 ( Dl-X-Y )*( D6-X ) +SKA* ( TDNS+2 .*X+2 . *Y )**2*! D6-X ) 
DFX4(X>Y)=(D3+Y)*6.*(D5+3. *X+2 . *Y )+SKB*( TDNS+2 . *X+2 . *Y ) 

1-2 . *SKB*( Dl-X-Y ) 

DFY4(X,Y)=4.*(05+3.*X+2.*Y)*(D3+Y)+(D5+3.*X+2.*Y)**2+SKB 
1*( TDNS+2. *X+2.*Y)-2.*SKB*( Dl-X-Y) 

F5(X,Y)=( (03+7.*Y-X)**7/10000.*(05+X+15.*Y)**15-SKA/10000. 
1*(TDNS+14.*Y)**14 
1*(D1-Y)*(D6-X-7.*Y)/10000. )/l.E10 
F6(X,Y)=(D4+X)*(D5+X+15.*Y)-SKB*(D6-X-7.*Y)*(D3+7.*Y-X) 
DFX5(X,Y)=(-(05+X+15.*Y)**15/10300.*7.*lD3+7.*Y-X)**6/100. 
1+(D3+7.*Y-X) **7/10000. *15. *(D5+X+15.*Y ) **14/1 00 . +SKA/1000000 . * 
41TDNS+ 

214.*Y)**14*(D1-Y) J/1.E08 

DFY5(X,Y)=(7.*(D3+7.*Y-X)**6/10000. *7.*(D5+X+15 .*Y )**15 
1+15./10000 .*!D5+X+15.*Y)**14 

1 *15.*(D3+7.*Y-X )**7+SKA/10000 . *( TDNS+14 . *Y )**14*( D6-X-7. *Y )-SKA - 

2/10000. *14. *( 

2TDNS + 14.*Y ) **13*14. *( Dl-Y )*( D6-X-7.*Y ) -SKA/10 000 . *( TDNS+14 . *Y )**14- 
4*( Dl-Y )*( -7. ) )/l .E10 

DFX6(X,Y)=(05+X+15.*Y)+(D4+X)+SKB*< ( D3+7. *Y-X ) + ( D6-X-7 . *Y ) ) 
DFY6(X,Y)=15.*(D4+X)-7.*SKB*( ( D6-X-7 . *Y )- C D3+7 . *Y-X ) ) 

C CAL. THE TOTAL AMOUNT 

IF ( ( IJ . EQ. 2 ) .OR . ( IJ.EQ.3 ) ) GO TO 51 

DO 5 IA=1,7 

DNS( IA )=DNS( IA 1/1000 . 

5 CONTINUE 
51 TONS-O. 

DO 1 IA=1,I 
TDNS=TCNS+DNS(IA) 

1 CCNTIN'JE 
Dl-DNSt 1 ) 

D2=DNS( 2 ) 

D3=DNS( 3) 

D4=DMS(4) 

D5=CNS(5) 
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0119480 


D6=DN5( 6 ) 

0119500 


D7=DHS(7) 

0119520 


SKA=SK1(1) 

0119540 


SKB=SK1 ( 2 ) 

0119560 


IF (IJ.EQ.l) GO TO 61 

0119580 


IF (IJ.EQ.2) GO TO 62 

0119600 


X : XY( 2 1/1000 . 

0119620 


Y=XY(1)/1000. 

0119640 


DO 2 IA=1,500 

0119660 


D=0FX5(X,Y)*DFY6fX,Y)-DFX6(X,Y)*DFY5(X,Y) 

0119680 


DX=(F6(X,Y)»DFY5(X,Y)-F5(X,Y)*DFY6CX,Y) )/D 

0119700 


0Y=(F5(X»Y)#DFX6(X»Y)-F6(X,Y )#DFX5( X» Y) )/D 

0119720 


X=X+DX 

0119740 


Y=Y+OY 

0119760 


AX=ABS(DX) 

0119780 


AY-ABSf DY ) 

0119300 


IF C (AX. LT. 0.001). AND. (AY. LT. 0.001)') GO TO 9 

0119820 

2 

CONTINUE 

0119840 


WHITE (6, 101) 

0119360 

101 

FCRMATt IX, 'SOLVE THE EQU. EQUATION FAIL AFTER 500 ITERATIONS’) 

0119680 


RETURN 

0119900 

9 

XY( 1 )-Y*1000 . 

0119920 


XY(2)=X#1000. 

0119940 

61 

X=XY(2) 

0119960 


Y=XY(1) 

0119930 


DO 3 IA=1,500 

0120000 


D=0FX1(X,Y)#DFY2(X,Y)-0FX2(X,Y)*DFY1(X,Y) 

0120020 


DX=(F2(X,Y)*DFYl<X,Y)--Fl(X,Y)*DFY2tX,Y) )/D 

0120040 


DY=(F1(X,Y)#DFX2(X,Y)-F2(X,Y)*DFX1(X,Y ) )/0 

0120060 


X=X+DX 

0120C80 


Y=Y+DY 

0120100 


AX=ABSfDX) 

0120120 


AY=ABS(DY) 

0120140 


IFUAX. LT. 0.001). AND. (AY. LT. 0.001)) GO TO 91 

0120160 

3 

CONTINUE 

0120180 


WRITE( 6 , 101 ) 

0120200 


RETURN 

0120220 

91 

XY(1)=Y 

0120240 


XY(2)=X 

0120260 


RETURN 

0120280 

62 

X=XY(1) 

0120300 


Y=XY(2) 

0120320 


DO 4 IA=1,500 

0120340 


D=DFX3(X,Y)*DFY4(X,Y)-DFX4(X,Y)#DFY3(X,Y) 

0120360 


DX=(F4(X,Y)*0FY3(X,Y)-F3(X,Y)*0FY4tX. Y) )/D 

0120380 


DY=(F3(X,Y)#DFX4(X,Y)-F4(X,Y)#DFX3(X,Y) )/D 

0120400 


X=X+DX 

0120420 


Y=Y+DY 

’0120440 


AX=ABS(DX) 

0120460 


AY=ABSfOY) 

0120480 


IF((AX. LT. 0.001). AND. (AY. LT. 0.001)) GO TO 92 

01205C0 

4 

CONTINUE 

0120520 


WRITE ( 6 1 101 ) 

0120540 


RETURN 
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0120560 92 XY(1)=X 

0120560 XY(2)=Y 

C120600 RETURN 

0120620 END 

0120640 SUBROUTINE VINEW ( M , V ,Z ,TK , POP , PPH2 , PF02 , PPH20 , PPCO , XO ) 

0120660 CCMMON/CATAL/SRO ,SA,CU,CL,ALFA,SN > FCONST > AREAF >DKC 

0120680 1 R-8 . 314 

0120700 ERR = 0 . 005 

0120720 CLA=CL 

0120740 EZ=1.261-.00025*TK 

0120760 SR=SRC*EXP( 3650 . *( 1 ./TK-1 ./450 . ) ) 

0120780 SIO=. 2327*( PP02*F0P)**. 6*1 PPH20*P0P)**.4377*EXP( -6652. /TK) 

012C600 C=SIO*SA*CU*CL 

0120320 EX=11.85*.0066*PPCO*POP*EXP(9190.*< 1./TK-1./450. ) ) 

0120340 A=ALC5( PPH2/PFH20*( PF02*P0P )**0 .5 ) 

0120860 C1=CLA*SA*CU*. 000053 

0120380 D=R*TK/SN/FCONST 

0120900 B=EZ+D*A 

C120920 DA=D/ALFA 

0120940 CDL=CNC/AREAF*(PP02*P0P) 

0120960 IF (M.EQ.2) GO TO 2 

01209S0 V=B-DA*ALCG( Z/C ) -Z*SR-EX*ALOG( Z/Cl )-D*ALOG( CDL/( CDL-Z ) ) 

0121000 GO TO 6 

0121020 2 Z=X0 

0121040 3 CONTINUE 

0121060 DO 5 1=1,50 

01210S0 FZ=Z*SR+DA*AL03( Z/C )+V-B+EX*ALGG( Z/Cl )+D*ALOG( CD L/( CDL-Z ) ) 

0121100 OFZ=SR+DA/Z+EX/Z+D/( CDL-Z) 

0121120 DZ=FZ/DFZ 

C121140 Z=Z-DZ 

C121160 4 IF (Z.LE.O.) GO TO 7 

0121180 IF (ABS(DZ).LT.ERR) GO TO 6 

01212C0 5 CONTINUE 

0121220 7 Z=X0 

0121240 DO 8 1=1,50 

0121260 GFZ=( DA*ALOG( Z/C )+V-B+EX*ALCG( Z/Cl )+D*ALOG(CDL/( CDL-Z )) )/( -SR ) 

0121260 GZ=Z 

0121300 Z=3FZ 

0121320 IF (Z.LE.O.) GO TO 19 

0121340 IF ( ABS( ( GZ-Z )/( Z+GZ ) ) . LT . ERR ) GO TO 6 

0121360 8 CONTINUE 

0121330 KRITEt 6 ,201 ) 

0121400 201 FORMAT! IX, 'CURRENT DENSITY LOOPING') 

0121420 GO TO 6 

0121440 19 ERR=ERR+0 . 001 
0121460 GO TO 2 

0121480 6 CONTINUE 

0121500 RETURN 

0121520 END 

0121540 SUBROUTINE KREF ( DNSR ,DNSF ,DX,DY ,DQ , PX , TCO ,THZ ,Z , POUT ,TCOUT 

0121560 1 .THOUT , SI ,DP1 , IFUEL ) 

0121580 REAL K0.MH.K1.K2 

0121600 COMMON FO ,F1 , F2 , F3 , F4 ,F5 ,F6 ,MH ,CG1 ,CG2 ,CG3,CG4,CG5,CG6 

0121620 COMMON/REP/ KO , EA ,RHOB , EPS.DZZ 
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0121640 

0121660 

0121660 

0121700 

0121720 

0121740 

0121760 

0121780 

0121800 

0121820 

0121840 

0121860 

0121880 

0121900 

0121920 

0121940 

0121960 

0121980 

0122000 

0122020 

0122040 

0122060 

0122080 

0122100 

0122120 

0122140 

0122160 

0122180 

0122200 

0122220 

0122240 

0122260 

0122230 

0122300 

0122320 

0122340 

0122360 

0122330 

0122400 

0122420 

0122440 

0122460 

0122480 

0122500 

0122520 

0122540 

0122560 

0122580 

'0122600 

0122620 

0122640 

0122660 

0122680 

0122700 


COMHON/ADDRE/ D1 ,D2 ,03,S,0P, P 
COMMON /KM/ KM 
COMMON /FCG/ F7.CG7 

DIMENSION Xl( 50),TAX2( 50),XE2( 50),TC( 50), 

1TH ( 50 ) >XE( 50 ) >TA(50 ) ,TAK1( 50),TAK2( 50), 

2XF ( 50 ) , XCOMPt 50,7), XMCOMPC 50,7),XCf 50,7), 

3XU1( 50 ) > TP( 50),FL(7),C(7),WM(7),COM(7), CGCOMPt 7 ) ,TCGC( 7 ) 
DIMENSION REN( 50) 

DIMENSION P( 50) 

DIMENSION 0NSR( 7) ,DNSF( 7) 

DATA EFROR/O . 01/ 

DATA XCOMP/ 350*0./ 

IDEBUS=0 

01=DX 

D2=DT 

D3=DQ 

DP=DP1 

S=S1 

P(1) = PX 

F1=DNSR(1) 

F2=DNSR(3) 

F3=DNSR(4) 

F4=DNSR( 6 ) 

F5-DNSR( 5 ) 

F6=DNSR( 7) 

F7-DNSR( 2 ) 

F0=F1+F2+F3+F4+F5+F6+F7 

CG1=DNSF(1) 

CG2=DNSF(3) 

CG3-DNSF( 4 ) 

CG4=0NSF(6) 

CG5=DNSF(5) 

CG6=DNSF( 7) 

CG7=DNSF(2) 

MH=CG1+CG2+CG3+CG4+CG5+CG6+CG7 
X1<1) = 0. 

XE2( 1 ) = 0 . 

TCO=( TCO-273 . 16 )*9. /5 . +32 . 

THZ=( THZ-Z73. 16 )*9./5. +32 . 

KIN3=0 

CO=(Fl*(l.-Xl(l))*P(l))/( . 730 2* ( F0 + 2.*X1(1)*F1)*( TCO+460 . ) ) 
U0=(4.*.7302*(F0+2.*X1(1)*F1 )*( TCO+46 0. ) )/( 3 . 1415927*P( 1 ) 
1*(D2**2-D1**2)*EPS) 

FIRST ASSUMPTION -- THO 
X=.9 

AX=Z*RH03*KO*P( 1 )*F1*( 2 .*F0+X*( 2.*F1-F0 ) )/( Z.*F0*( FO + 2 .*X*F1 )* 
1U0*C0) 

T=- . 9045*EA/AIOG( X/AX ) -460 . 

TK2=K2(T) 

TX2=X2(X,TK2) 

TDH1 = -DH1( T ) 

TDH2=-DH2(T) 

AT=Fl*TDHl*X+< F3+X*F1 )*TDH2*TX2 
FL(1) = F1*(1.-X) 



0122720 

0122740 

0122760 

0122780 

0122800 

0122820 

0122840 

0122660 

0122880 

0122900 

0122920 

0122940 

0122960 

0122930 

0123000 

0123020 

0123040 

0123060 

0123080 

0123100 

0123120 

0123140 

0123160 

0123180 

0123200 

0123220 

0123240 

0123260 

0123280 

0123300 

0123320 

0123340 

0123360 

0123380 

0123400 

0123420 

0123440 

0123460 

0123480 

0123500 

0123520 

0123540 

0123560 

0123580 

0123600 

0123620 

0123640 

0123660 

'0123680 

0123700 

0123720 

0123740 

0123760 

0123780 


FL( 2 )=F2-TX2«< F3+X»F1 ) 

F L ( 3)=F3+X*F1+TX2*( F3+X *F1) 
FL(4)=F4-2.»X#F1-TX2«(F3+X «F1) 
FL15)=F5+4.*X*F1+TX2*(F3+X *F1) 

FLt 6 ) = F6 
FL(7) = F7 
TF=F( X ) 

CALL CMASStC.FL.TF) 

TVIS-VIS( CiT i IFUEL ) 

TTHC=THCt C i T . IFUEL 1 
THI=HI(TVIS,TTHC) 

CGCOMPt 1 )=CG1 
CGCOMPt 2 )=CG2 
CGCOMPt 3 )=CG3 
CGCOMPt 4 )=CG4 
CGCOMPt 5 )=CG5 
CGCOMPt 6 )=CG6 
CGCOMPt 7) =CG7 
CALL CMASSt TCGC iCGCOMP >MH ) 

CGVIS-VISt TCGCiTHZ i IFUEL ) 

CGTHC=THC( TCGC, THZ, IFUEL) 

CGHTCP=HTCP( TCGC > THZ ) 

THO-HOt CGVIS > CGTHC »CGHTCP i Z iTHZ >RE ) 
TUI=UI(THI,THO,T) 

AZ=3 . 1415927*TUI*Z*D2/2 . 

CALL COMPINt COM ( X1 1 1 ) , XE2 ( 1 ) ) 

TFCP=FCPtT .COM ) 

AW=( MH#CGHTCP#TFCP )/AZ+MH#CGHTCP-TFCP 
THO-( AW*THZ+2 . #TFCP#TCO+AY )/( AW+2 .»TFCP) 
IF(THO.LT.TCO) THO=TCO+250. 

IDEA=1 
IHOPE=l 
75 1 = 1 

TAKlt 1 ) = 0 . 

TAK2(1) = 0. 

XUltl)=0. 

TAt 1 )=TCO 
TAX2t 1 ) = 0 . 

XF(l)=F0+2.*Fl*Xltl) 

DO 80 J=1 , 7 
XCOMPt 1 ( J )=COM( J ) 

XHCOMPt 1 , J )=XCOMP( 1 , J )/XF< 1 ) 

80 CONTINUE 
UMM1=0. 

DO 81 J=1 , 7 

WMMl=WMMl + < XCOMPt 1 , J )*WM( J ) )/XF( 1 ) 

81 CONTINUE 

DO 82 J=1 ( 7 

XCt 1, J )=( XCOMPt l,J)*WMt J))/tXF<l)*WMMl) 

82 CONTINUE 
TCt 1 )=TCO 
THt 1 )=THO 

SECOND ASSUMPTION — TCtI+1) 

72 TCt 1 + 1 )=TC( I ) 



0123800 

0123820 

0123340 

0123860 

0123880 

0123900 

0123920 

0123940 

0123960 

0123930 

0124000 

0124020 

0124040 

0124060 

0124080 

0124100 

0124120 

0124140 

0124160 

0124180 

0124200 

0124220 

0124240 

0124260 

0124280 

0124300 

0124320 

0124340 

0124360 

0124380 

0124400 

0124420 

0124440 

0124460 

0124480 

0124500 

0124520 

0124540 

0124560 

0124530 

0124600 

0124620 

0124640 

0124660 

0124680 

0124700 

0124720 

0124740 

0124760 

0124780 

0124800 

0124820 

0124840 

0124860 


70 TA(I+1) = (TC(I )+TC( 1 + 1 ) )/2 . 

TAK1(I+1)=K1(TA(I+1) ) 

TAK2(I+1)=K2(TA(I+1)) 

81=DZZ»RHOB*KO*P< I )/( U0*C0 ) 

B2=EXP( - . 9045+EA/I TA( 1 + 1 1 + 460. 1 ) 

Xl(I+l)=Xl(I)+Bl*E2*XnCOMP<I,l) 

XE2(I + 1)=X2(X1(I + 1),TAK2(I + 1)1 
XDH1 = -DH1(TA(I+D) 

XDH2=-DH2(TA(I+11) 

AL1=F1*XDH1*(X1(I+1)-X1( ID 
AL2=(F3+X1(I+1)*F1)*XDH2*(XEZ<I+1)-XE2(I) ) 

AL=AL1+AL2 

FL(1)=F1#(1.-X1(I+1)) 

FL(2)-F2-XE2(I+1)#(F3+X1(I+1)*F1) 

FL(3J=F3+X1( 1 + 1 )*F1+XE2( 1+1 )*( F3+XK 1 + 1 )*F1 ) 
FL(4)=F4-2.*X1(I+1)*F1-XE2(I+1)*(F3+X1(I+1)*F1) 

FL( 5 )=F5+4. »X1 (I+1)#F1+XE2(I+1)*( F3+XM 1+1 )»F1 ) 

FL( 6 )=F6 
FL( 7 )=F7 

XF(I + 1) = F(X1(I+1)) 

CALL CMASS(C,FL,XF( 1+1 ) ) 

DO 10 J = l,7 
XCOMPt I+1,J)=FL(J) 

XMCOHP(I+l,J)=XCOMP(I+l,J)/XF(I+l) 

XC(I+1,J)=C(J) 

10 CONTINUE 

XVIS=VIS( C,TA( 1+1 ) , IFUEL) 

XTHC=THC(C,TA(I+1), IFUEL) 

V=;:i1( 1 )*F1+U'M( 3 )*F2+UM(4 )*F3+Wf1(6 )«F4+KM( 5 )*F5+WM( 7 )*F6 +WM( 2 )*F7 
G!1V=( V*4. )/( 3. 1415927»( D2**2-D1**2 ) ) 

AHN=V/F0 

RHO=(AKW#P(I) )/< 0.7302#(T+460. ) ) 

DELP=( 1 . -EPS )/EPS#*3*GMV**2/DP/4 . 18E+08/RHO»( 150 .*( 1 . -EPS )*XVIS/DP 
1/GMV+l . 75 )*DZZ/2116 . 8 
Pt 1 + 1 ) = P( I l-DELP 
XHI=HI( XVIS,XTHC ) 

XCGVIS-VIS( TCGC , TH( 1 1 , IFUEL 1 
XCGTHC=THC( TCGC ,TH( 1 1 , IFUEL) 

XGHTCP=HTCP( TCGC, TH( 111 
XHO=HO(XCGVIS,XCGTHC,XGHTCP,Z,TH< I), RE) 

REN( 1+1 )=RE 

XU1(I+1)=UI(XHI,XH0,TA(I+1 ) ) 

AM-( 3. 1415927#XU1( 1+1 )*0ZZ*D2 )/2. 

XFCP=FCP(TA(I),FL) 

AN= ( KH*XGHTCP*XFCP )/AM-XFCP+MH*XGHTCP 

TH( 1+1 )=TH( I )*( AN+2 . »XFCP )/AN-2.*XFCP*TC( I l/AN-AL/AN 

TP( 1+1 ) =TH ( 1+1 )*( AM-MH*XGHTCP)/Af1+TH( I >*( AM+MH»XGHTCP )/AM-TC( I ) 

C TEST SECOND ASSUMPTION 

EE=ABS< ( TP( 1 + 1 )-TC( 1 + 1 ) )/TC( 1+1 ) ) 

IF ( IDEBUG . ME . 0 ) 

1WRITE( 6 , 2020 ) TP( 1 + 1 ) ,TC( 1 + 1 ) , EE 
2020 FORMATt 1 OTP= 1 , 1PE15 . 7,5X, 1 TC= 1 , E15. 7,5X, 1 EE= 1 , E15. 7) 

IFfEE.LE. ERROR ) GO TO 71 

IF( KING. LE . 15 ) TC( 1 + 1 ) = ( TC( 1 + 1 )+TP( 1 + 1 ) )/2 . 



0124380 

0124900 

0124920 

0124940 

0124960 

0124930 

0125000 

0125020 

0125040 

0125060 

0125060 

0125100 

0125120 

0125140 

0125160 

0125180 

0125200 

0125220 

0125240 

0125260 

0125280 

0125300 

0125320 

0125340 

0125360 

0125380 

0125400 

0125420 

0125440 

0125460 

0125480 

0125500 

0125520 

0125540 

0125560 

0125530 

0125600 

0125620 

0125640 

0125660 

0125680 

0125700 

0125720 

0125740 

0125760 

0125780 

0125800 

0125820 

’0125640 

0125360 

0125880 

0125900 

0125920 

0125940 


TC2=TC(I+1) 

EETC2=TH( 1 + 1 )*( AM-MH*XGHTCP )/AM+TH( I )*( AM+MH#XGHTCP )/AM- - 

1TC( I )-TC2 

IF( KING . LE . 15 ) GO TO 97 

IF( EETC1 .NE . EETC2 ) TC3=TC2-< EETC2/< EETC2-EETC1 ) )*ITC2-TC1 ) 

IF( EETC1 . EQ.EETC2 ) TC3=( TC2+TP( 1 + 1 ) 1/2 
TC( 1 + 1 )=TC3 
97 KIN3=KIN3+1 
TC1=TC2 
EETC1=EETC2 

IF ( KING . GE . 40 ) GO TO 959 
GO TO 70 
71 CONTINUE 

TC( 1 + 1 )=TP( 1+1 ) 

AA=I 

AAA=AA*DZZ 

IF(AAA.GE.Z) GO TO 73 
1=1 + 1 
KIHS=0 
GO TO 72 

73 N=I+1 

TEST FIRST ASSUMPTION 

AB=ABS( ( TH( 1+1 l-THZl/THZ) 

IF(AS.LE. 0.001) GO TO 74 

TH02=TH0 

THZ2=TH( N ) 

IF ( IDEA . LT . 2 ) TH03=TH0+THZ-THZ2 

IF( IDEA -EQ . 2 ) TH03= ( TH01-TH02 )/( THZ1-THZ2 )*l THZ-THZ2 )+TH02 

IDEA=2 

THZ1=TH( N ) 

TH01=TH0 

TH0=TH03 

IFITHO. LT.TCO) IHOPE=IHOPE+l 
IF ( THO . LT . TOO ) TH0=TC0+5Q. 

IF ( IH3PE.EQ.5) GO TO 975 
GO TO 75 

74 CONTINUE 
IFIK.EQ.1) L=N 
GO TO 954 

975 CONTINUE 
URITEl 6 > 976 ) THO 

976 FCRMATC 1H1 > ' »»* INSUFFICIENT COMB. GAS HEAT CAPACITY ’ / 
l'OTHO=' ,F17.3, 'THIS IS LESS THAN TCO'/'CRAISE THZ AND/OR COM3. 
2GAS FLOW RATES' ) 

GO TO 954 

959 NRITE( 6 > 958 ) KING, I 

958 FORMATUHl,' LOOPING ON TC FOR ' ,14 ITERATIONS IN INCRM',14) 
954 CONTINUE 

TCO=( TCO-32 . 1*5. /9. +273. 16 
THZ=( THZ-32 . )»5 ./9. +273. 16 
FCUT=P( N ) 

TCCUT=(TC(N)-32. )»5./9. +273. 16 
THOUT=(TH(l)-32. )*5./9. +273. 16 
DNSRl 1 )=XCOMPt N> 1 ) 
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0125960 DNSR(3)=XCOMP(N,2) 

0125980 DNSR ( 4 ) =XCOMP ( N > 3 ) 

0126000 DNSR(5)=XC0MP(N,5) 

0126020 DNSR ( 6 )=XCOMP( N,4 ) 

0126040 0NSR(7)=XC0NP(N,6) 

0126060 RETURN 

01260S0 END 

0163100 SUBROUTINE CMASSt C, FL,F ) 

0163200 DIMENSION C< 7 ) ,WM( 7 ) , FL< 7 I 

0163300 COMMON /U’M/ WM 

0163400 KMM z (FL(l)*WM(l) + FL(2)*U71<3) + FL(3)*Wm4)+WM(6)*FL(4) + FL(5)*WM(5) + - 

0163500 lFL(6)*nT(7) + FL(7)*HM(2) )/F 

0163600 Cl 1 ) = FL( 1 )*KM< 1 )/( F#WMM ) 

0163700 Cl 2 ) = FL( 2 )*UMt 3 )/( F#WMM ) 

0163800 Cl 3 )=FL( 3 )*U'M( 4 )/( F*WMM ) 

0163900 C(4) z FL(4)#KM(6)/( F*WMM ) 

0164000 C(5) = FL(5)»WM(5)/1F*KMM) 

0164100 Cl 6 ) = FL( 6 )#WM( 7 )/( F*WMM ) 

0164200 Cl 7) = FL( 7 )#HM( 2 )/( F#NMM ) 

0164300 RETURN 

0164400 END 

0164500 SUBROUTINE COMPINl COM, XI ,X2 ) 

0164600 REAL MH 

0164700 COMMON FO , FI , F2 , F3 , F4 , F5 , F6 ,MH ,CG1 ,CG2 , CG3 ,CG4 ,CG5 , CG6 

0164S00 COMMON/ADORE/ D1 ,02 ,03, S ,DP, P 

0164900 COMMON /FCG/F7.CG7 

0165000 DIMENSION CCM17) 

0165100 DO 8 J=l,7 

0165200 IF(J.EQ.l) COM! J ) = F1*( 1 . -XI ) 

0165300 IF1J.EQ.2) COM! J ) Z F2-X2*< F3+X1*F1 ) 

0165400 IF1J.EQ.3) COM( J ) Z F3+X1*F1+X2*( F3+X1#F1 ) 

0165500 IF1J.EQ.4) COMU ) z F4-2 . #X1#F1-X2*( F3+X1#F1 ) 

0165600 IF1J.EQ.5) COM( J ) z F5+4 . #X1*F1+X2*( F3+X1*F1 ) 

0165700 IF1J.EQ.6) COM(J)=F6 

0165800 IF1J.EQ.7) C0M(J)=F7 

0165900 8 CONTINUE 

0166000 RETURN 

0166100 END 

0166200 FUNCTION DH11T) 

0166300 DHl=-2.7285E-03*(T#*2)+12.698»T+7.002E+04 

0166400 RETURN 

0166500 END 

0166600 FUNCTION 0H2(T) 

01667C0 DH2=2 . 328Q#T-18111 .4 

0166800 RETURN 

0166900 END 

0167000 FUNCTION F(X1) 

0167100 REAL MH 

0167200 COMMON FO , FI , F2 , F3 , F4 , F5 , F6 ,MH , CGI ,CG2 ,CG3 , CG4 ,CG5 ,CG6 

0167300 COMMGN/ADDRE/ D1 ,02 ,D3 , S , DP , P 

0167400 COMMON /FCG/F7.CG7 

0167500 F=F0+2.*X1#F1 

0167600 RETURN 

0167700 END 



0167800 

0167900 

0168000 

0168100 

0168200 

0168300 

0166400 

0168500 

0168600 

0168700 

0168600 

0168900 

0169000 

0169100 

0169200 

0169300 

0169400 

0169500 

0169600 

0169700 

0169800 

0169900 

0170000 

C1701C0 

0170200 

0170300 

0170400 

0170500 

0170600 

0170700 

0170S00 

0170900 

0171000 

0171100 

0171200 

0171300 

0171400 

0171500 

0171600 

0171700 

0171800 

0171900 

0172000 

0172100 

0172200 

0172300 

0172400 

0172500 

'0172600 

0172700 

0172SC0 

0172900 

0173000 

0173100 


FUNCTION FCP(T,COM) 

DIMENSION COM( 7) ,A(4,7) 

DATA A/5. 34, 6. 39E-03,0. ,0. ,6. 60,6.675-04,0. ,0. ,10. 34,1. 52E-03.0. , 
1-6 . 33420E + 05.8 . 22 , 8. 3E-05.4 . 136E-07 , 0 . ,6.62 .4.5E-04, 0 . , 0 . ,6.5, 
25.56E-04.0. ,0. ,6 . 732 ,8 . 36E-03 ,5 .53E-09, 0 . 0/ 

TP=T +460 . 

FCP-0 . 

DO 9 1=1,7 

FCP=FCP+COM(I>*< A(1,I)+A(2,I)*TP+A(3,I)*TP**2+A(4,I)/(TP**2>) 

9 CONTINUE 
RETURN 
END 

FUNCTION HI(VIS.THC) 

REAL MH 

DIMENSION KM( 7 ) 

COMMON F0,F1,F2,F3,F4,F5,F6,(1H,CG1 >CG2 >CG3,CG4,CG5 >CG6 

COMMCN/ADDR E/ D1 ,D2 ,03 , S ,DP , P 

COMMON/KM/ KM 

COMMON /FCG/F7.CG7 

HI = 0 . 

V=(KM(1)*F1+WM< 3)*F2+WM(4)*F3+WM(6)*F4+WM(5)*F5+WM(7)*F6 
1+KM( 2 )*F7) 

GMV= ( V*4 . )/(3.1415927*(D2**2.-Dl**2. )) 

HID=( .813*(GMV*0P/VIS)**.9)*EXP(-6.*DP/(DZ-D1) ) 

HI=( HI0*THC )/( D2-D1 ) 

RETURN 

END 

FUNCTION HO(VIS,THC,HTCP,Z,T,RE) 

REAL MH 

DIMENSION KM( 7 ) 

COMMON F0,F1,F2,F3,F4,F5,F6,MH , CGI , CG2 ,CG3 ,CG4 ,CG5 , CG6 
COMMCN/ADDR E/ D1 ,D2 ,D3,S ,DP, P 
COMMON /FCG/F7.CG7 
COMMON/KM/ KM 

AMN=( CG1*KM C 1 ) +CG2*KM( 3 ) +CG3*KM( 4 ) +CG4*KM( 6 ) +CG5*WM ( 5 ) +CG6*KM( 7 ) 
1+CG7+WM( 2 ) )/MH 
H0=0 . 

G=MH/( S**2-( 3. 1415927*03**2 )/4. ) 

DE=4.*(S**2-3. 1415927*03**2-/4. )/( 3 . 1415927*03+4 . *S ) 
RE=(DE*G*AMW)/VIS 
PR=( HTCP*VIS )/( THC*AMW 1 
RHO=(AMU*P)/(0.7302*(T+460. ) ) 

GR = ( Z**3 )*( RH0**2 1*4 . 18E08*100 . /( VIS**2 1 
IFIRE. GE. 10000 . ) GO TO 300 
IFIRE. LE. 2100. ) GO TO 200 
2100<RE<10000 

H02 1= ( 1 . *THC/DE )*f 2100 . ** . 45 )*SQRT( PR )*( DE/Z )** . 4*( S/D 3 )** . 8 
1*GR**.05 

H010 = ( . 02*THC/DE )*(10000.**.8)*(PR**.333)*( S/D 3 )**.53 
SLO?E = (HO10-HO21 1/(10000.-2100. 1 
H0=K021 + SL0FE*( RE-2100. 1 
RETURN 

RE<=2100 

200 HO J ( 1 . 02*THC/DE )*( RE**.45 1*SQRT( PR )*( DE/Z )**.4*(S/D3 )**.8 



128 


0173200 

0173300 

0173400 C 

0173500 

0173600 

0173700 

0173720C 

0173740C 

C173760C 

0173780C 

0173800C 

0173S20C 

0173640 

0173860 

0173880 

0173900 

0173920 

0173940 

0173960 

0173930 

0174000 

0174020C 

C174040C 

0174060C 

01 74030C 

0174100C 

0174120C 

0175400 

0175500 

0175600 

0175700 

0175600 

0175900 

0176000 

0176100 

0176200 

0176300 

0176400 

0176500 

0176600 

0176700 

0176800 

0176900 

0177000 

0177100 

0177200 

0177300 

0177400 

0177500 

0177600 

0177700 

0177800 

0177900 

0173000 


1*GR**.05 

RETURN 

RE>=10000 

300 HO=( .02*THC/DE)*(RE**.8)*< PR**. 333 )*( S/D 3 )**.53 
RETURN 
END 

FUNCTION HTCP(CM.T) 

DIMENSION Cm7),C(7),A(4,7),WN(7) 

DATA A/5.34,6.39E-03,0. ,0. ,6 .60,6 .67E-04 , 0 . ,0. , 10 . 34, 1 . 52E-03, 0 
1-6.3342E+05.8.22, 8. 3E-05.4. 136E-07, 0 . ,6 .62 .4.5E-04, 0 . , 0 . ,6 . 5, 
25 .56E-04 . 0 . > 0 . ,6 . 732 , 8. 36E-03.5.53E-09, 0 ./ 

DATA WM/16 . i 28 . ,44. ,18. ,2. ,28. ,32./ 

TC-CM( 1 )/W1( 1 )+CMC 2 )/UM< 2 J+CM( 3 )/NM( 3 )+CM( 4 )/WM( 4 ) +CM( 5 )/WMt 5 ) 
1+CM( 6 )/WM( 6 )+CM( 7 )/Uf1( 7 ) 

C( 1 )=CM( 1 )/W“( 1 J/TC 
C(2)=CM(2)/WM(2)/TC 
C( 3 ) =CMC 3 )/l.’M( 3 )/TC 
C( 4 ) = CM ( 4 )/MM( 4 )/TC 
CC5)=CM(5)/kM(5)/TC 
C(6)=CM(6)/W1(6)/TC 
C( 7 J =CMC 7 )/KM( 7 )/TC 
TP=T+460. 

HTCP=0. 

DO 1 1=1,7 

1 HTCP=HTCP+C(I)*(A(1,I)+A(2,I J*TP+A( 3,I)*TP**2+A(4,I )/TP**2 ) 
RETURN 
END 

REAL FUNCTION K1(T) 

TP=T+460. 

B=-8.7153Ef08/(TP**3)+5.2409E+06/(TP**2. J-4.6299E+04/TP+27.849 
K1=EXP( B ) 

RETURN 

END 

REAL FUNCTION K2(T) 

TP=T+460. 

B=-9.0283E + 08/( TP**3 )+3. 5603E+06/( TP**2 )+4. 366 2 E+ 03/TP- 3. 0526 
K2 = EXP( B ) 

RETURN 

END 

FUNCTION THC( C,T , IJ ) 

DIMENSION C( 7 ) , A( 2 , 7 ) 

CCMMON/THCC/ A 
THC=0. 

DO 5 1=1,7 

THC=THC+C(I)*(A(1,I)*T+A(2,m 
5 CONTINUE 
RETURN 
END 

FUNCTION UI(HI,HO,T) 

REAL MH 

COMNON/AODRE/ D1 ,D2,03,S,DP,P 
COMMON /FCG/F7.CG7 
DLM=(03-D2 l/ALOG! 03/D2 J 
R=0 . 005 



0178100 

0178200 

0178300 

0178400 

0178500 

01786C0 

0178700 

017S800 

0178900 

0179000 

0179100 

0179200 

0179300 

C179400 

0179500 

0179600 

0179700 

0179800 

0179900 

0180000 

0180100 

0180200 

0180300 

016C400 

0130500 

0180600 

0160700 

0180800 

0180900 

0181000 

0181100 

0181200 

0181300 

01814C0 

0181500 

0181600 

0181700 

0181800 

0181900 

01S2000 

0182100 

0182200 

0182300 

0182400 

0182500 

0182600 

0182700 

0132800 

'0182900 


THMET=4.659E-03*T+6.248 

UI=l./(l./HI+D2/(D3*HO)+< (D3-D2 )«D2 )/(THMET*DU1)+R ) 

RETURN 

END 

FUNCTION VIS( C iT, IJ ) 

DIMENSION At 2 ,7 ) iC( 7 ) 

COMMCN/VIFC/ A 
VIS=0. 

5 DO 4 1=1,7 

VIS=VI5+C( I )*l At 1 ,1 )«T+At 2,1)) 

4 CONTINUE 
RETURN 
END 

FUNCTION X2t X ,K2 ) 

REAL NH 
REAL K2 

COMMON F0,F1,F2,F3,F4,F5,F6,MH,CG1,CG2,CG3,CG4,CG5,CG6 
COM.ilON/ADDRE/ 01 , D2 ,D3 , S ,DP , P 
COMMON /FCG/F7.CG7 
A=( K2-1 . )*( F3+X*F1 )#*2 

B=t F3+X*F1)*(2.*X*F1*K2-K2*F2-K2*F4-5.«X*F1-F3-F5) 
C=K2*F2*F4-2.*Fl*=2*K2*X-(F3+X*Fl)*tF5+4.*Fl»X) 

X2 = ( -5-SORT ( B** 2 -4 . *A*C ) )/( 2 .*A ) 

IF(X2.LT.-1. .OR.X2.GT.O. 1 X2=( -B+SQRTt B**2-4.*A*C ) )/( 2 . *A ) 

RETURN 

END 

SOPFC T0PFC=443. ,UT=0.8,CD=325. , 

SEND 

SINIT DNSM=1216. ,0. ,1.360,21.8,166. ,0 . ,0 . ,TAT=293. , PAT=1 . ,SMRA=3. ,POPR=5.0 
SEND 

SCCNDT IFUEL=1,ERR=0.01,IP=2,I=7,EXT=100. ,WAT=0 . 015 , EXA=100 . , 

SEND 

SREPEN ZH=6. ,DX1=0. ,0X2=0 . 15 ,DX3=0 . 1667, K0=1.040E+04 , EA=20000 . ,RHOB=80. 

,EPS=0.487,S=0.25,DP=0.00328,DZZ=0.25, 

SEND 

SHEATX CN=1.3,U=48825.1,HA(7)=0.2,HA(10)=0.2, 

SEND 

SHEFDC NFH=2 , NRH=5 , BSPAC=1 . ,ODTH= . 0625 , PITCH= . 0833 ,CLH= . 0208 , IDSH= . 833 , 
IDTH=. 04667, FLOAR=. 001716, SURFC=. 1466, CLENH=2. ,S1TS2=0. 5, DTH=0. 7, 

SEND 

SFDSHH DPD=1 . 18, 0 . , AHRN=0 .66 , . 0 , APPD=69. ,0. ,CLEPD=5. 91 ,0 . , 

NTPD=1 , 0 , 

SEND 

SFDFUH NTAF=140 , FULE =1. 42, WIDAF=. 009744, NPFU=3365 , 

NTAA=40 , AIRL=1 . ,UIOAA= . 00515 
SEND 

SCATAI SRO= .44 ,SA=400 . , CU= . 15 ,CL= . 75 , ALFA= .5, SN=2 . , FCONST= 96500 . , 
DKC=2.4E5, 

SEND 
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